%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{The sadists package} %\VignetteKeyword{distributions} %\VignettePackage{sadists} \documentclass[10pt,a4paper,english]{article} % front matter%FOLDUP \usepackage[hyphens]{url} \usepackage{amsmath} \usepackage{amsfonts} % for therefore \usepackage{amssymb} % for theorems? \usepackage{amsthm} \theoremstyle{plain} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{corollary}[theorem]{Corollary} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \newtheorem{conjecture}[theorem]{Conjecture} \newtheorem{example}{Example}[section] \theoremstyle{remark} \newtheorem*{remark}{Remark} \newtheorem*{caution}{Caution} \newtheorem*{note}{Note} % see http://tex.stackexchange.com/a/3034/2530 \PassOptionsToPackage{hyphens}{url}\usepackage{hyperref} \usepackage{hyperref} \usepackage[square,numbers]{natbib} %\usepackage[authoryear]{natbib} %\usepackage[iso]{datetime} %\usepackage{datetime} %compactitem and such: \usepackage[newitem,newenum,increaseonly]{paralist} \makeatletter \makeatother %\input{sr_defs.tex} \usepackage{sadists} %\providecommand{\sideWarning}[1][0.5]{\marginpar{\hfill\includegraphics[width=#1\marginparwidth]{warning}}} % see: https://stat.ethz.ch/pipermail/r-help/2007-November/144810.html \newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}} \let\proglang=\textsf \let\code=\texttt \newcommand{\CRANpkg}[1]{\href{https://cran.r-project.org/package=#1}{\pkg{#1}}} \newcommand{\sadists}{\CRANpkg{sadists}\xspace} \newcommand{\PDQutils}{\CRANpkg{PDQutils}\xspace} % knitr setup%FOLDUP <<'preamble', include=FALSE, warning=FALSE, message=FALSE, cache=FALSE>>= library(knitr) # set the knitr options ... for everyone! # if you unset this, then vignette build bonks. oh, joy. #opts_knit$set(progress=TRUE) opts_knit$set(eval.after='fig.cap') # for a package vignette, you do want to echo. # opts_chunk$set(echo=FALSE,warning=FALSE,message=FALSE) opts_chunk$set(warning=FALSE,message=FALSE) #opts_chunk$set(results="asis") opts_chunk$set(cache=FALSE,cache.path="cache/sadists") #opts_chunk$set(fig.path="figure/",dev=c("pdf","cairo_ps")) opts_chunk$set(fig.path="figure/sadists",dev=c("png","pdf")) #opts_chunk$set(fig.width=4.0,fig.height=6,dpi=200) # see http://stackoverflow.com/q/23419782/164611 opts_chunk$set(fig.width=4.0,fig.height=2.5,out.width="5in",out.height="3in",dpi=144) # doing this means that png files are made of figures; # the savings is small, and it looks like shit: #opts_chunk$set(fig.path="figure/",dev=c("png","pdf","cairo_ps")) #opts_chunk$set(fig.width=4,fig.height=4) # for figures? this is sweave-specific? #opts_knit$set(eps=TRUE) # this would be for figures: #opts_chunk$set(out.width='.8\\textwidth') # for text wrapping: options(width=64,digits=2) #opts_chunk$set(size="tiny") opts_chunk$set(size="small") opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=36,blank=TRUE)) # from the environment # compiler flags! library(sadists) @ %UNFOLD % SYMPY preamble%FOLDUP %\usepackage{graphicx} % Used to insert images %\usepackage{adjustbox} % Used to constrain images to a maximum size \usepackage{color} % Allow colors to be defined \usepackage{enumerate} % Needed for markdown enumerations to work %\usepackage{geometry} % Used to adjust the document margins \usepackage{amsmath} % Equations \usepackage{amssymb} % Equations %\usepackage[utf8]{inputenc} % Allow utf-8 characters in the tex document %\usepackage[mathletters]{ucs} % Extended unicode (utf-8) support \usepackage{fancyvrb} % verbatim replacement that allows latex %\usepackage{grffile} % extends the file name processing of package graphics % to support a larger range % The hyperref package gives us a pdf with properly built % internal navigation ('pdf bookmarks' for the table of contents, % internal cross-reference links, web links for URLs, etc.) \usepackage{hyperref} %\usepackage{longtable} % longtable support required by pandoc >1.10 \definecolor{orange}{cmyk}{0,0.4,0.8,0.2} \definecolor{darkorange}{rgb}{.71,0.21,0.01} \definecolor{darkgreen}{rgb}{.12,.54,.11} \definecolor{myteal}{rgb}{.26, .44, .56} \definecolor{gray}{gray}{0.45} \definecolor{lightgray}{gray}{.95} \definecolor{mediumgray}{gray}{.8} \definecolor{inputbackground}{rgb}{.95, .95, .85} \definecolor{outputbackground}{rgb}{.95, .95, .95} \definecolor{traceback}{rgb}{1, .95, .95} % ansi colors \definecolor{red}{rgb}{.6,0,0} \definecolor{green}{rgb}{0,.65,0} \definecolor{brown}{rgb}{0.6,0.6,0} \definecolor{blue}{rgb}{0,.145,.698} \definecolor{purple}{rgb}{.698,.145,.698} \definecolor{cyan}{rgb}{0,.698,.698} \definecolor{lightgray}{gray}{0.5} % bright ansi colors \definecolor{darkgray}{gray}{0.25} \definecolor{lightred}{rgb}{1.0,0.39,0.28} \definecolor{lightgreen}{rgb}{0.48,0.99,0.0} \definecolor{lightblue}{rgb}{0.53,0.81,0.92} \definecolor{lightpurple}{rgb}{0.87,0.63,0.87} \definecolor{lightcyan}{rgb}{0.5,1.0,0.83} % commands and environments needed by pandoc snippets % extracted from the output of `pandoc -s` \DefineShortVerb[commandchars=\\\{\}]{\|} \DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}} % Add ',fontsize=\small' for more characters per line \newenvironment{Shaded}{}{} \newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{{#1}}}} \newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{{#1}}} \newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}} \newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}} \newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}} \newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}} \newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}} \newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{{#1}}}} \newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{{#1}}} \newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}} \newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{{#1}}} \newcommand{\RegionMarkerTok}[1]{{#1}} \newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}} \newcommand{\NormalTok}[1]{{#1}} % Define a nice break command that doesn't care if a line doesn't already % exist. \def\br{\hspace*{\fill} \\* } % Math Jax compatability definitions \def\gt{>} \def\lt{<} % Pygments definitions \makeatletter \def\PY@reset{\let\PY@it=\relax \let\PY@bf=\relax% \let\PY@ul=\relax \let\PY@tc=\relax% \let\PY@bc=\relax \let\PY@ff=\relax} \def\PY@tok#1{\csname PY@tok@#1\endcsname} \def\PY@toks#1+{\ifx\relax#1\empty\else% \PY@tok{#1}\expandafter\PY@toks\fi} \def\PY@do#1{\PY@bc{\PY@tc{\PY@ul{% \PY@it{\PY@bf{\PY@ff{#1}}}}}}} \def\PY#1#2{\PY@reset\PY@toks#1+\relax+\PY@do{#2}} \expandafter\def\csname PY@tok@gd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.00,0.00}{##1}}} \expandafter\def\csname PY@tok@gu\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.50,0.00,0.50}{##1}}} \expandafter\def\csname PY@tok@gt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.27,0.87}{##1}}} \expandafter\def\csname PY@tok@gs\endcsname{\let\PY@bf=\textbf} \expandafter\def\csname PY@tok@gr\endcsname{\def\PY@tc##1{\textcolor[rgb]{1.00,0.00,0.00}{##1}}} \expandafter\def\csname PY@tok@cm\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}} \expandafter\def\csname PY@tok@vg\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}} \expandafter\def\csname PY@tok@m\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@mh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@go\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.53,0.53}{##1}}} \expandafter\def\csname PY@tok@ge\endcsname{\let\PY@it=\textit} \expandafter\def\csname PY@tok@vc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}} \expandafter\def\csname PY@tok@il\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@cs\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}} \expandafter\def\csname PY@tok@cp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.74,0.48,0.00}{##1}}} \expandafter\def\csname PY@tok@gi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.63,0.00}{##1}}} \expandafter\def\csname PY@tok@gh\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}} \expandafter\def\csname PY@tok@ni\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.60,0.60,0.60}{##1}}} \expandafter\def\csname PY@tok@nl\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.63,0.00}{##1}}} \expandafter\def\csname PY@tok@nn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}} \expandafter\def\csname PY@tok@no\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.00,0.00}{##1}}} \expandafter\def\csname PY@tok@na\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.49,0.56,0.16}{##1}}} \expandafter\def\csname PY@tok@nb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@nc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}} \expandafter\def\csname PY@tok@nd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}} \expandafter\def\csname PY@tok@ne\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.82,0.25,0.23}{##1}}} \expandafter\def\csname PY@tok@nf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}} \expandafter\def\csname PY@tok@si\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}} \expandafter\def\csname PY@tok@s2\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@vi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}} \expandafter\def\csname PY@tok@nt\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@nv\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}} \expandafter\def\csname PY@tok@s1\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@sh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@sc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@sx\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@bp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@c1\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}} \expandafter\def\csname PY@tok@kc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@c\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}} \expandafter\def\csname PY@tok@mf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@err\endcsname{\def\PY@bc##1{\setlength{\fboxsep}{0pt}\fcolorbox[rgb]{1.00,0.00,0.00}{1,1,1}{\strut ##1}}} \expandafter\def\csname PY@tok@kd\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@ss\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}} \expandafter\def\csname PY@tok@sr\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}} \expandafter\def\csname PY@tok@mo\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@kn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@mi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@gp\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}} \expandafter\def\csname PY@tok@o\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}} \expandafter\def\csname PY@tok@kr\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@s\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@kp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@w\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.73,0.73}{##1}}} \expandafter\def\csname PY@tok@kt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.69,0.00,0.25}{##1}}} \expandafter\def\csname PY@tok@ow\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}} \expandafter\def\csname PY@tok@sb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \expandafter\def\csname PY@tok@k\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}} \expandafter\def\csname PY@tok@se\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.13}{##1}}} \expandafter\def\csname PY@tok@sd\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}} \def\PYZbs{\char`\\} \def\PYZus{\char`\_} \def\PYZob{\char`\{} \def\PYZcb{\char`\}} \def\PYZca{\char`\^} \def\PYZam{\char`\&} \def\PYZlt{\char`\<} \def\PYZgt{\char`\>} \def\PYZsh{\char`\#} \def\PYZpc{\char`\%} \def\PYZdl{\char`\$} \def\PYZhy{\char`\-} \def\PYZsq{\char`\'} \def\PYZdq{\char`\"} \def\PYZti{\char`\~} % for compatibility with earlier versions \def\PYZat{@} \def\PYZlb{[} \def\PYZrb{]} \makeatother % Exact colors from NB \definecolor{incolor}{rgb}{0.0, 0.0, 0.5} \definecolor{outcolor}{rgb}{0.545, 0.0, 0.0} % Prevent overflowing lines due to hard-to-break entities \sloppy % Setup hyperref package \hypersetup{ breaklinks=true, % so long urls are correctly broken across lines colorlinks=true, urlcolor=blue, linkcolor=darkorange, citecolor=darkgreen, } % Slightly bigger margins than the latex defaults %\geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in} %UNFOLD %UNFOLD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % document incantations%FOLDUP \begin{document} \title{The \pkg{sadists} package} \author{Steven E. Pav}% \thanks{\email{shabbychef@gmail.com}}} %\date{\today, \currenttime} \maketitle %UNFOLD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{abstract}%FOLDUP The \sadists package includes `dpqr' functions for some obscure distributions, mostly involving sums and ratios of (non-central) chi-squares, chis, and normals. \end{abstract}%UNFOLD \section{Introduction}%FOLDUP The \sadists package provides density, distribution, quantile and random generation functions (the `dpqr' functions) for some obscure distributions. For all of these, the `dpq' functions are approximated via the Edgeworth and Cornish-Fisher expansions. As such, this package is a showcase for the capabilities of the \PDQutils package, which does the heavy lifting once the cumulants have been computed. \cite{PDQutils-Manual} It should be noted that the functions provided by \sadists do \emph{not} recycle their distribution parameters against the \Rcode{x, p, q} or \Rcode{n} parameters. This is in contrast to the common \Rlang idiom, and may cause some confusion. This is mostly for reasons of performance, but also because some of the distributions have vector-valued parameters; recycling over these would require the user to provide \emph{lists} of parameters, which would be unpleasant. First, a function which will evaluate the `dpq' functions versus random draws of the variable: <<'setup',echo=TRUE,eval=TRUE>>= require(ggplot2) require(grid) DEFAULT_NOBS <- 2^17 testf <- function(dpqr,nobs,...) { rv <- sort(dpqr$r(nobs,...)) data <- data.frame(draws=rv,pvals=dpqr$p(rv,...)) text.size <- 6 # sigh eat <- c(0.005,0.01,0.05,0.10,0.20,0.35) eq <- sort(unique(c(eat,0.5,1-eat))) sumdata <- data.frame(thrpv=eq,emppv=quantile(data$pvals,eq), thrrv=dpqr$q(eq,...),emprv=quantile(data$draws,eq)) # http://stackoverflow.com/a/5688125/164611 p1 <- ggplot(data, aes(x=draws)) + geom_line(aes(y = ..density.., colour = 'Empirical'), stat = 'density') + stat_function(fun = function(x) { dpqr$d(x,...) }, aes(colour = 'Theoretical')) + geom_histogram(aes(y = ..density..), alpha = 0.3) + scale_colour_manual(name = 'Density', values = c('red', 'blue')) + theme(text=element_text(size=text.size)) + labs(title="Density (tests dfunc)") # Q-Q plot p2 <- ggplot(sumdata, aes(x=thrrv, y=emprv)) + geom_point() + geom_abline(slope=1,intercept=0,colour='red') + theme(text=element_text(size=text.size)) + labs(title="Q-Q plot (tests qfunc)",x="Theoretical",y="Sample", caption="Data subsampled to selected quantiles.") # empirical CDF of the p-values; should be uniform p3 <- ggplot(sumdata, aes(x=thrpv, y=emppv)) + geom_point() + geom_abline(slope=1,intercept=0,colour='red') + theme(text=element_text(size=text.size)) + labs(title="P-P plot (tests pfunc)",x="Theoretical",y="Sample", caption="Data subsampled to selected quantiles.") # Define grid layout to locate plots and print each graph pushViewport(viewport(layout = grid.layout(2, 2))) print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1:2)) print(p2, vp = viewport(layout.pos.row = 2, layout.pos.col = 1)) print(p3, vp = viewport(layout.pos.row = 2, layout.pos.col = 2)) } @ %UNFOLD \section{Sum of (non-central) chi-squares to a power}%FOLDUP Let $X_i \sim \chi^2_{\nu_i}\left(\delta_i\right)$ be independent non-central chi-square variates, where $\delta_i$ are the non-centrality parameters and $\nu_i$ are the degrees of freedom. Let $w_i, p_i$ be given constants. Then $$ Y = \sum_i w_i X_i^{p_i} $$ follows a weighted sum of non-central chi-squares to a power distribution. This is not a common distribution. However, its cumulants can be easily computed, so the `pdq' functions can be approximated by classical expansions. Moreover, its CDF and quantile functions can be used to compute those of the doubly non-central F, and it is related to the upsilon distribution. <<'sumchisqpow',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the sum of chi-squares to a power distribution.">>= require(sadists) wts <- c(-1,1,3,-3) df <- c(100,200,100,50) ncp <- c(0,1,0.5,2) pow <- c(1,0.5,2,1.5) testf(list(d=dsumchisqpow,p=psumchisqpow,q=qsumchisqpow,r=rsumchisqpow),nobs=DEFAULT_NOBS,wts,df,ncp,pow) @ %UNFOLD \section{K-prime distribution}%FOLDUP \label{sec:kprime_distribution} Let $X_i \sim \chi^2_{\nu_i}$ be chi-square random variables with $\nu_i$ degrees of freedom for $i=1,2$, independent of $Z\sim\mathcal{N}\left(0,1\right)$, a standard normal. Suppose $a, b$ are given constants. Then $$ Y = \frac{b Z + a \sqrt{X_1 / \nu_1}}{\sqrt{X_2 / \nu_2}} $$ follows a K-prime distribution with degrees of freedom $\asrowvec{\nu_1, \nu_2}$ and parameters $a, b$. \cite{PoitLecoutre2010} Depending on these four parameters, the K-prime generalizes the following: \begin{itemize} \item The normal distribution, when $b=1, a=0, \nu_2 = \infty$. \item The Lambda-prime distribution (see \secref{lamprime_distribution}), when $b=1, a\ne 0, \nu_2 = \infty$. \item The (central) t-distribution, when $b=1, a=0, \nu_2 < \infty$. \item The square-root of the F-distribution, when $b = 0,a = 1$. \item The (central) chi-distribution, when $b=0, a = 1, \nu_2 = \infty$. \end{itemize} <<'kprime',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the K-prime distribution.">>= require(sadists) v1 <- 50 v2 <- 80 a <- 0.5 b <- 1.5 testf(list(d=dkprime,p=pkprime,q=qkprime,r=rkprime),nobs=DEFAULT_NOBS,v1,v2,a,b) @ %UNFOLD \section{Lambda prime distribution}%FOLDUP \label{sec:lamprime_distribution} Let $X \sim \chi^2_{\nu}$ be a chi-square random variable with $\nu$ degrees of freedom, independent of $Z\sim\mathcal{N}\left(0,1\right)$, a standard normal. Then $$ Y = Z + t \sqrt{X/\nu} $$ follows a Lambda-prime distribution with parameter $t$ and degrees of freedom $\nu$. \cite{Lecoutre1999} It is a special case of the K-prime distribution (\secref{kprime_distribution}) and of the upsilon distribution (\secref{upsilon_distribution}). <<'lambdap',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the Lambda-prime distribution.">>= require(sadists) df <- 50 ts <- 1.5 testf(list(d=dlambdap,p=plambdap,q=qlambdap,r=rlambdap),nobs=DEFAULT_NOBS,df,ts) @ %UNFOLD \section{Upsilon distribution}%FOLDUP \label{sec:upsilon_distribution} Let $X_i \sim \chi^2_{\nu_i}$ be independent central chi-square variates, where $\nu_i$ are the degrees of freedom. Let $Z\sim\mathcal{N}\left(0,1\right)$ be a standard normal, independent of the $X_i$. Let $t_i$ be given constants. Then $$ Y = Z + \sum_i t_i \sqrt{X_i} $$ follows an upsilon distribution with parameter $\asrowvec{t_1,t_2,\ldots,t_k}$ and degrees of freedom $\asrowvec{\nu_1,\nu_2,\ldots,\nu_k}$. <<'upsilon',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the upsilon distribution.">>= require(sadists) df <- c(30,50,100,20,10) ts <- c(-3,2,5,-4,1) testf(list(d=dupsilon,p=pupsilon,q=qupsilon,r=rupsilon),nobs=DEFAULT_NOBS,df,ts) @ %UNFOLD \section{Doubly non-central F distribution}%FOLDUP \label{sec:dnf} The doubly non-central F distribution generalizes the F distribution to the case where the denominator chi-square is non-central. For $i=1,2$, let $X_i \sim \chi^2_{\nu_i}\left(\delta_i\right)$ be independent non-central chi-square variates, where $\delta_i$ are the non-centrality parameters and $\nu_i$ are the degrees of freedom. Then $$ Y = \frac{X_1/\nu_1}{X_2/\nu_2} $$ follows a doubly non-central F distribution. <<'dnf',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the doubly non-central F distribution.">>= require(sadists) df1 <- 40 df2 <- 80 ncp1 <- 1.5 ncp2 <- 2.5 testf(list(d=ddnf,p=pdnf,q=qdnf,r=rdnf),nobs=DEFAULT_NOBS,df1,df2,ncp1,ncp2) @ %UNFOLD \section{Doubly non-central t distribution}%FOLDUP \label{sec:dnt} The doubly non-central t distribution generalizes the t distribution to the case where the denominator chi-square is non-central. Let $X_2 \sim \chi^2_{\nu_2}\left(\delta_2\right)$ be a non-central chi-square variate, with non-centrality parameter $\delta_2$ and $\nu_2$ degrees of freedom. Let $X_2$ be independent of $Z$, a standard normal. Then $$ Y = \frac{Z + \delta_1}{\sqrt{X_2/\nu_2}} $$ follows a doubly non-central t distribution with degrees of freedom $\nu_2$ and non-centrality parameters $\delta_1, \delta_2$. The square of a doubly non-central t is, up to scaling, a doubly non-central F, see \secref{dnf}. <<'dnt',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the doubly non-central t distribution.">>= require(sadists) df <- 75 ncp1 <- 2 ncp2 <- 3 testf(list(d=ddnt,p=pdnt,q=qdnt,r=rdnt),nobs=DEFAULT_NOBS,df,ncp1,ncp2) @ %UNFOLD \section{Doubly non-central Beta distribution}%FOLDUP \label{sec:dnbeta} The doubly non-central Beta distribution generalizes the Beta distribution to the case where the denominator chi-square is non-central. For $i=1,2$, let $X_i \sim \chi^2_{\nu_i}\left(\delta_i\right)$ be independent non-central chi-square variates, where $\delta_i$ are the non-centrality parameters and $\nu_i$ are the degrees of freedom. Then $$ Y = \frac{X_1}{X_1 + X_2} $$ follows a doubly non-central Beta distribution. Note that $$ F = \frac{\nu_2}{\nu_1}\frac{Y}{1-Y} $$ follows a doubly non-central F distribution. The `PDQ' functions use this relationship. <<'dnbeta',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the doubly non-central Beta distribution.">>= require(sadists) df1 <- 40 df2 <- 80 ncp1 <- 1.5 ncp2 <- 2.5 testf(list(d=ddnbeta,p=pdnbeta,q=qdnbeta,r=rdnbeta),nobs=DEFAULT_NOBS,df1,df2,ncp1,ncp2) @ %UNFOLD \section{Doubly non-central Eta distribution}%FOLDUP \label{sec:dneta} The doubly non-central Eta distribution is to the doubly non-central Beta what the doubly non-central t is to the doubly non-central F. Let $X_2 \sim \chi^2_{\nu_2}\left(\delta_2\right)$ be a non-central chi-square variate, with non-centrality parameter $\delta_2$ and $\nu_2$ degrees of freedom. Let $X_2$ be independent of $Z$, which is normal with mean $\delta_1$ and unit standard deviation. Then $$ Y = \frac{Z}{\sqrt{Z^2 + X_2}} $$ follows a doubly non-central Eta distribution with degrees of freedom $\nu_2$ and non-centrality parameters $\delta_1, \delta_2$. The square of a doubly non-central Eta is a doubly non-central Beta, see \secref{dnbeta}. Note that $$ t = {\sqrt{\nu_2}}\frac{Y}{\sqrt{1-Y^2}} $$ follows a doubly non-central t distribution. The `PDQ' functions use this relationship. <<'dneta',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the doubly non-central Eta distribution.">>= require(sadists) df <- 100 ncp1 <- 0.5 ncp2 <- 2.5 testf(list(d=ddneta,p=pdneta,q=qdneta,r=rdneta),nobs=DEFAULT_NOBS,df,ncp1,ncp2) @ %UNFOLD \section{Sum of logs of (non-central) chi-squares}%FOLDUP \label{sec:sumlogchisq} Let $X_i \sim \chi^2_{\nu_i}\left(\delta_i\right)$ be independent non-central chi-square variates, where $\delta_i$ are the non-centrality parameters and $\nu_i$ are the degrees of freedom. Let $w_i$ be given constants. Then $$ Y = \sum_i w_i \log{X_i} $$ follows a weighted sum of log of non-central chi-squares distribution. This is not a common distribution. However, its cumulants can easily be computed. \cite{pav2015lnc} <<'sumlogchisq',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the sum of log of chi-squares distribution.">>= require(sadists) wts <- c(5,-4,10,-15) df <- c(100,200,100,50) ncp <- c(0,1,0.5,2) testf(list(d=dsumlogchisq,p=psumlogchisq,q=qsumlogchisq,r=rsumlogchisq),nobs=DEFAULT_NOBS,wts,df,ncp) @ %UNFOLD \section{Product of (non-central) chi-squares to a power}%FOLDUP Let $X_i \sim \chi^2_{\nu_i}\left(\delta_i\right)$ be independent non-central chi-square variates, where $\delta_i$ are the non-centrality parameters and $\nu_i$ are the degrees of freedom. Let $p_i$ be given constants. Then $$ Y = \prod_i X_i^{p_i} $$ follows a product of non-central chi-squares to a power distribution. This is not a common distribution. The `PDQ' functions are computed using a transform of the logs of chi-squares distribution, see \secref{sumlogchisq}. <<'prodchisqpow',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the product of chi-squares to a power distribution.">>= require(sadists) df <- c(100,200,100,50) ncp <- c(0,1,0.5,2) pow <- c(1,0.5,2,1.5) testf(list(d=dprodchisqpow,p=pprodchisqpow,q=qprodchisqpow,r=rprodchisqpow),nobs=DEFAULT_NOBS,df,ncp,pow) @ %UNFOLD \section{Product of doubly non-central F variates}%FOLDUP Let $X_j \sim F_{\nu_{1,j},\nu_{2,j}}\left(\delta_{1,j},\delta_{2,j}\right)$ be independent doubly non-central F variates, where $\delta_{i,j}$ are the non-centrality parameters and $\nu_{i,j}$ are the degrees of freedom. Then $$ Y = \prod_j X_j $$ follows a product of doubly non-central Fs distribution. This is not a common distribution. The `PDQ' functions are computed using a transform of the logs of chi-squares distribution, see \secref{sumlogchisq}. <<'proddnf',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the product of doubly non-central Fs distribution.">>= require(sadists) df1 <- c(10,20,5) df2 <- c(1000,500,150) ncp1 <- c(1,0,2.5) ncp2 <- c(0,1.5,5) testf(list(d=dproddnf,p=pproddnf,q=qproddnf,r=rproddnf),nobs=DEFAULT_NOBS,df1,df2,ncp1,ncp2) @ %UNFOLD \section{Product of normal variates}%FOLDUP Let $Z_j \sim \mathcal{N}\left(\mu_{j},\sigma_{j}^2\right)$ be independent normal variates with means $\mu_{j}$ and variances $\sigma^2_{j}$. Then $$ Y = \prod_j X_j $$ follows a product of normals distribution. This is not a common distribution. When the coefficients of variation, $\sigma_j / \mu_j$ are large for some of the variates, the approximations given in this package tend to break down. <<'prodnormal',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the product of normals distribution.">>= require(sadists) mu <- c(100,20,5) sigma <- c(10,2,0.2) testf(list(d=dprodnormal,p=pprodnormal,q=qprodnormal,r=rprodnormal),nobs=DEFAULT_NOBS,mu,sigma) @ %UNFOLD %\section{Sum of generalized gamma variates}%FOLDUP %Let $X_i \sim G\left(a_i, d_i, p_i\right)$ %be independently distributed generalized gamma variates with parameters %$a_i, d_i$ and inverse power $p_i$. \cite{stacey1962} %Let $w_i$ be given constants. Then %$$ %Y = \sum_i w_i X_i %$$ %follows a sum of generalized gammas distribution. %This is not a common distribution. %<<'sumgengamma',echo=TRUE,eval=TRUE,fig.cap="Confirming the dpqr functions of the sum of generalized gammas distribution.">>= %require(sadists) %wts <- c(8,-10) %a <- c(10,2) %d <- c(2,1) %pow <- c(0.5,1.5) %testf(list(d=dsumgengamma,p=psumgengamma,q=qsumgengamma,r=rsumgengamma),nobs=DEFAULT_NOBS,wts,a,d,pow) %@ %%UNFOLD \nocite{paolella2007intermediate,Lecoutre1999,PoitLecoutre2010} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % bibliography%FOLDUP \nocite{markowitz1952portfolio,markowitz1999early,markowitz2012foundations} %\bibliographystyle{jss} %\bibliographystyle{siam} %\bibliographystyle{ieeetr} \bibliographystyle{plainnat} %\bibliographystyle{acm} \bibliography{sadists,rauto} %UNFOLD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %\appendix%FOLDUP %\section{Confirming the scalar Gaussian case} %\begin{example}%FOLDUP %To sanity check %\theoremref{theta_asym_var_gaussian}, %consider the $\nlatf = 1$ Gaussian case. In this case, %$$ %\fvech{\pvsm} = \asvec{1,\pmu,\psig^2 + \pmu^2},\quad\mbox{and}\quad %\fvech{\minv{\pvsm}} = %\asvec{1 + \frac{\pmu^2}{\psig^2},- \frac{\pmu}{\psig^2},\oneby{\psig^2}}. %$$ %Let $\smu, \ssig^2$ be the unbiased sample estimates. By well %known results \cite{spiegel2007schaum}, $\smu$ and $\ssig^2$ are independent, %and have asymptotic variances of $\psig^2/\ssiz$ and $2\psig^4/\ssiz$ %respectively. By the delta method, the asymptotic variance of %$\Unun \fvech{\svsm}$ and $\fvech{\minv{\svsm}}$ can be computed as %\begin{equation} %\begin{split} %\VAR{\Unun \fvech{\svsm}} &\rightsquigarrow %\oneby{\ssiz}\qform{\twobytwossym{\psig^2}{0}{2\psig^4}}{\twobytwo{1}{2\pmu}{0}{1}},\\ %&= \oneby{\ssiz}\twobytwossym{\psig^2}{2\pmu\psig^2}{4\pmu^2\psig^2 + %2\psig^4}. %\label{eqn:gauss_confirm_theta} %\end{split} %\end{equation} %\begin{equation} %\begin{split} %\VAR{\fvech{\minv{\svsm}}} &\rightsquigarrow %\oneby{\ssiz}\qform{\twobytwossym{\psig^2}{0}{2\psig^4}}{\twobythree{\frac{2\pmu}{\psig^2}}{-\frac{1}{\psig^2}}{0}{-\frac{\pmu^2}{\psig^4}}{\frac{\pmu}{\psig^4}}{-\frac{1}{\psig^4}}},\\ %&= %\oneby{\ssiz}\gram{\twobythree{2\psnr}{-\frac{1}{\psig}}{0}{-\sqrt{2}\psnr^2}{\sqrt{2}\frac{\psnr}{\psig}}{-\frac{\sqrt{2}}{\psig^2}}},\\ %&= \oneby{\ssiz}\threebythreessym{2\psnr^2\wrapParens{2 + \psnr^2}}{-\frac{2\psnr}{\psig}\wrapParens{1+\psnr^2}}{2\frac{\psnr^2}{\psig^2}}{\frac{1 + %2\psnr^2}{\psig^2}}{-\frac{2\psnr}{\psig^3}}{\frac{2}{\psig^4}}. %\label{eqn:gauss_confirm_itheta} %\end{split} %\end{equation} %Now it remains to compute $\VAR{\Unun \fvech{\svsm}}$ via %\theoremref{theta_asym_var_gaussian}, and then %\VAR{\fvech{\minv{\svsm}}} via %\theoremref{inv_distribution}, and confirm they match the values above. %This is a rather tedious computation best left to a computer. Below is %an excerpt of an iPython notebook using Sympy \cite{PER-GRA:2007,sympy} %which performs this computation. %This notebook is available online. \cite{SEP2013example} %% SYMPY from here out%FOLDUP %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}1}]:} \PY{c}{\PYZsh{} confirm the asymptotic distribution of Theta} %\PY{c}{\PYZsh{} for scalar Gaussian case.} %\PY{k+kn}{from} \PY{n+nn}{\PYZus{}\PYZus{}future\PYZus{}\PYZus{}} \PY{k+kn}{import} \PY{n}{division} %\PY{k+kn}{from} \PY{n+nn}{sympy} \PY{k+kn}{import} \PY{o}{*} %\PY{k+kn}{from} \PY{n+nn}{sympy.physics.quantum} \PY{k+kn}{import} \PY{n}{TensorProduct} %\PY{n}{init\PYZus{}printing}\PY{p}{(}\PY{n}{use\PYZus{}unicode}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{,} \PY{n}{wrap\PYZus{}line}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{,} \PYZbs{} %\PY{n}{no\PYZus{}global}\PY{o}{=}\PY{n+nb+bp}{True}\PY{p}{)} %\PY{n}{mu} \PY{o}{=} \PY{n}{symbols}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{\PYZbs{}}\PY{l+s}{mu}\PY{l+s}{\PYZsq{}}\PY{p}{)} %\PY{n}{sg} \PY{o}{=} \PY{n}{symbols}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{\PYZbs{}}\PY{l+s}{sigma}\PY{l+s}{\PYZsq{}}\PY{p}{)} %\PY{c}{\PYZsh{} the elimination, duplication and U\PYZus{}\PYZob{}\PYZhy{}1\PYZcb{} matrices:} %\PY{n}{Elim} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{3}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %\PY{n}{Dupp} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{4}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %\PY{n}{Unun} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{2}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %\PY{k}{def} \PY{n+nf}{Qform}\PY{p}{(}\PY{n}{A}\PY{p}{,}\PY{n}{x}\PY{p}{)}\PY{p}{:} %\PY{l+s+sd}{\PYZdq{}\PYZdq{}\PYZdq{}compute the quadratic form x\PYZsq{}Ax\PYZdq{}\PYZdq{}\PYZdq{}} %\PY{k}{return} \PY{n}{x}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)} \PY{o}{*} \PY{n}{A} \PY{o}{*} \PY{n}{x} %\end{Verbatim} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}2}]:} \PY{n}{Theta} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{2}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{n}{mu}\PY{p}{,}\PY{n}{mu}\PY{p}{,}\PY{n}{mu}\PY{o}{*}\PY{o}{*}\PY{l+m+mi}{2} \PY{o}{+} \PY{n}{sg}\PY{o}{*}\PY{o}{*}\PY{l+m+mi}{2}\PY{p}{]}\PY{p}{)} %\PY{n}{Theta} %\end{Verbatim} %\texttt{\color{outcolor}Out[{\color{outcolor}2}]:} %\begin{equation*} %\left[\begin{matrix}1 & \mu\\\mu & \mu^{2} + \sigma^{2}\end{matrix}\right] %\end{equation*} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}3}]:} \PY{c}{\PYZsh{} compute tensor products and } %\PY{c}{\PYZsh{} the derivative d vech(Theta\PYZca{}\PYZhy{}1) / d vech(Theta)} %\PY{c}{\PYZsh{} see also \theoremref{inv_distribution}} %\PY{n}{Theta\PYZus{}Theta} \PY{o}{=} \PY{n}{TensorProduct}\PY{p}{(}\PY{n}{Theta}\PY{p}{,}\PY{n}{Theta}\PY{p}{)} %\PY{n}{iTheta\PYZus{}iTheta} \PY{o}{=} \PY{n}{TensorProduct}\PY{p}{(}\PY{n}{Theta}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)}\PY{p}{,}\PY{n}{Theta}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)}\PY{p}{)} %\PY{n}{theta\PYZus{}i\PYZus{}deriv} \PY{o}{=} \PY{n}{Elim} \PY{o}{*} \PY{p}{(}\PY{n}{iTheta\PYZus{}iTheta}\PY{p}{)} \PY{o}{*} \PY{n}{Dupp} %\end{Verbatim} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}4}]:} \PY{c}{\PYZsh{} towards \theoremref{theta_asym_var_gaussian}} %\PY{n}{DTTD} \PY{o}{=} \PY{n}{Qform}\PY{p}{(}\PY{n}{Theta\PYZus{}Theta}\PY{p}{,}\PY{n}{Dupp}\PY{p}{)} %\PY{n}{D\PYZus{}DTTD\PYZus{}D} \PY{o}{=} \PY{n}{Qform}\PY{p}{(}\PY{n}{DTTD}\PY{p}{,}\PY{n}{theta\PYZus{}i\PYZus{}deriv}\PY{p}{)} %\PY{n}{iOmega} \PY{o}{=} \PY{n}{Qform}\PY{p}{(}\PY{n}{D\PYZus{}DTTD\PYZus{}D}\PY{p}{,}\PY{n}{Unun}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)}\PY{p}{)} %\PY{n}{Omega} \PY{o}{=} \PY{l+m+mi}{2} \PY{o}{*} \PY{n}{iOmega}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)} %\PY{n}{simplify}\PY{p}{(}\PY{n}{Omega}\PY{p}{)} %\end{Verbatim} %\texttt{\color{outcolor}Out[{\color{outcolor}4}]:} %\begin{equation*} %\left[\begin{matrix}\sigma^{2} & 2 \mu \sigma^{2}\\2 \mu \sigma^{2} & 2 \sigma^{2} \left(2 \mu^{2} + \sigma^{2}\right)\end{matrix}\right] %\end{equation*} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}5}]:} \PY{c}{\PYZsh{} this matches the computation in \eqnref{gauss_confirm_theta}} %\PY{c}{\PYZsh{} on to the inverse:} %\PY{c}{\PYZsh{} actually use \theoremref{inv_distribution}} %\PY{n}{theta\PYZus{}i\PYZus{}deriv\PYZus{}t} \PY{o}{=} \PY{n}{theta\PYZus{}i\PYZus{}deriv}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)} %\PY{n}{theta\PYZus{}inv\PYZus{}var} \PY{o}{=} \PY{n}{Qform}\PY{p}{(}\PY{n}{Qform}\PY{p}{(}\PY{n}{Omega}\PY{p}{,}\PY{n}{Unun}\PY{p}{)}\PY{p}{,}\PY{n}{theta\PYZus{}i\PYZus{}deriv\PYZus{}t}\PY{p}{)} %\PY{n}{simplify}\PY{p}{(}\PY{n}{theta\PYZus{}inv\PYZus{}var}\PY{p}{)} %\end{Verbatim} %\texttt{\color{outcolor}Out[{\color{outcolor}5}]:} %\begin{equation*} %\left[\begin{matrix}\frac{2 \mu^{2}}{\sigma^{4}} \left(\mu^{2} + 2 \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{2 \mu^{2}}{\sigma^{4}}\\- \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{1}{\sigma^{4}} \left(2 \mu^{2} + \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}}\\\frac{2 \mu^{2}}{\sigma^{4}} & - \frac{2 \mu}{\sigma^{4}} & \frac{2}{\sigma^{4}}\end{matrix}\right] %\end{equation*} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}6}]:} \PY{c}{\PYZsh{} this matches the computation in \eqnref{gauss_confirm_itheta}} %\PY{c}{\PYZsh{} now check \conjectureref{theta_asym_var_gaussian}} %\PY{n}{conjec} \PY{o}{=} \PY{n}{Qform}\PY{p}{(}\PY{n}{Theta\PYZus{}Theta}\PY{p}{,}\PY{n}{Dupp}\PY{p}{)} %\PY{n}{e1} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{3}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{]}\PY{p}{)} %\PY{n}{convar} \PY{o}{=} \PY{l+m+mi}{2} \PY{o}{*} \PY{p}{(}\PY{n}{conjec}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)} \PY{o}{\PYZhy{}} \PY{n}{e1} \PY{o}{*} \PY{n}{e1}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)}\PY{p}{)} %\PY{n}{simplify}\PY{p}{(}\PY{n}{convar}\PY{p}{)} %\end{Verbatim} %\texttt{\color{outcolor}Out[{\color{outcolor}6}]:} %\begin{equation*} %\left[\begin{matrix}\frac{2 \mu^{2}}{\sigma^{4}} \left(\mu^{2} + 2 \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{2 \mu^{2}}{\sigma^{4}}\\- \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{1}{\sigma^{4}} \left(2 \mu^{2} + \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}}\\\frac{2 \mu^{2}}{\sigma^{4}} & - \frac{2 \mu}{\sigma^{4}} & \frac{2}{\sigma^{4}}\end{matrix}\right] %\end{equation*} %\begin{Verbatim}[commandchars=\\\{\}] %{\color{incolor}In [{\color{incolor}7}]:} \PY{c}{\PYZsh{} are they the same?} %\PY{n}{simplify}\PY{p}{(}\PY{n}{theta\PYZus{}inv\PYZus{}var} \PY{o}{\PYZhy{}} \PY{n}{convar}\PY{p}{)} %\end{Verbatim} %\texttt{\color{outcolor}Out[{\color{outcolor}7}]:} %\begin{equation*} %\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] %\end{equation*} %%UNFOLD %%% OLD%FOLDUP %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}1}]:} \PY{k+kn}{from} \PY{n+nn}{\PYZus{}\PYZus{}future\PYZus{}\PYZus{}} \PY{k+kn}{import} \PY{n}{division} %%\PY{k+kn}{from} \PY{n+nn}{sympy} \PY{k+kn}{import} \PY{o}{*} %%\PY{k+kn}{from} \PY{n+nn}{sympy.physics.quantum} \PY{k+kn}{import} \PY{n}{TensorProduct} %%\PY{n}{init\PYZus{}printing}\PY{p}{(}\PY{n}{use\PYZus{}unicode}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{,} \PY{n}{wrap\PYZus{}line}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{,} \PY{n}{no\PYZus{}global}\PY{o}{=}\PY{n+nb+bp}{True}\PY{p}{)} %%\PY{n}{mu} \PY{o}{=} \PY{n}{symbols}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{\PYZbs{}}\PY{l+s}{mu}\PY{l+s}{\PYZsq{}}\PY{p}{)} %%\PY{n}{sg} \PY{o}{=} \PY{n}{symbols}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{\PYZbs{}}\PY{l+s}{sigma}\PY{l+s}{\PYZsq{}}\PY{p}{)} %%\PY{c}{\PYZsh{} the elimination, duplication and U\PYZus{}\PYZob{}\PYZhy{}1\PYZcb{} matrices:} %%\PY{n}{Elim} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{3}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %%\PY{n}{Dupp} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{4}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %%\PY{n}{Unun} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{2}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{)} %%\PY{k}{def} \PY{n+nf}{quad\PYZus{}form}\PY{p}{(}\PY{n}{A}\PY{p}{,}\PY{n}{x}\PY{p}{)}\PY{p}{:} %%\PY{l+s+sd}{\PYZdq{}\PYZdq{}\PYZdq{}compute the quadratic form x\PYZsq{}Ax\PYZdq{}\PYZdq{}\PYZdq{}} %%\PY{k}{return} \PY{n}{x}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)} \PY{o}{*} \PY{n}{A} \PY{o}{*} \PY{n}{x} %%\end{Verbatim} %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}2}]:} \PY{n}{Theta} \PY{o}{=} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{2}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{n}{mu}\PY{p}{,}\PY{n}{mu}\PY{p}{,}\PY{n}{mu}\PY{o}{*}\PY{o}{*}\PY{l+m+mi}{2} \PY{o}{+} \PY{n}{sg}\PY{o}{*}\PY{o}{*}\PY{l+m+mi}{2}\PY{p}{]}\PY{p}{)} %%\PY{n}{Theta} %%\end{Verbatim} %%\texttt{\color{outcolor}Out[{\color{outcolor}2}]:} %%\begin{equation*} %%\left[\begin{matrix}1 & \mu\\\mu & \mu^{2} + \sigma^{2}\end{matrix}\right] %%\end{equation*} %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}3}]:} \PY{c}{\PYZsh{} compute tensor products and } %%\PY{c}{\PYZsh{} the derivative d vech(Theta\PYZca{}\PYZhy{}1) / d vech(Theta)} %%\PY{n}{Theta\PYZus{}Theta} \PY{o}{=} \PY{n}{TensorProduct}\PY{p}{(}\PY{n}{Theta}\PY{p}{,}\PY{n}{Theta}\PY{p}{)} %%\PY{n}{iTheta\PYZus{}iTheta} \PY{o}{=} \PY{n}{TensorProduct}\PY{p}{(}\PY{n}{Theta}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)}\PY{p}{,}\PY{n}{Theta}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)}\PY{p}{)} %%\PY{n}{theta\PYZus{}i\PYZus{}deriv} \PY{o}{=} \PY{n}{Elim} \PY{o}{*} \PY{p}{(}\PY{n}{iTheta\PYZus{}iTheta}\PY{p}{)} \PY{o}{*} \PY{n}{Dupp} %%\end{Verbatim} %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}4}]:} \PY{c}{\PYZsh{} towards thm \ref{theorem:theta_asym_var_gaussian}} %%\PY{n}{DTTD} \PY{o}{=} \PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{Theta\PYZus{}Theta}\PY{p}{,}\PY{n}{Dupp}\PY{p}{)} %%\PY{n}{D\PYZus{}DTTD\PYZus{}D} \PY{o}{=} \PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{DTTD}\PY{p}{,}\PY{n}{theta\PYZus{}i\PYZus{}deriv}\PY{p}{)} %%\PY{n}{iOmega} \PY{o}{=} \PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{D\PYZus{}DTTD\PYZus{}D}\PY{p}{,}\PY{n}{Unun}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)}\PY{p}{)} %%\PY{n}{Omega} \PY{o}{=} \PY{l+m+mi}{2} \PY{o}{*} \PY{n}{iOmega}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)} %%\PY{n}{simplify}\PY{p}{(}\PY{n}{Omega}\PY{p}{)} %%\end{Verbatim} %%\texttt{\color{outcolor}Out[{\color{outcolor}4}]:} %%\begin{equation*} %%\left[\begin{matrix}\sigma^{2} & 2 \mu \sigma^{2}\\2 \mu \sigma^{2} & 2 \sigma^{2} \left(2 \mu^{2} + \sigma^{2}\right)\end{matrix}\right] %%\end{equation*} %%This matches the computation in \eqnref{gauss_confirm_theta}. On to the %%inverse: %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}5}]:} \PY{c}{\PYZsh{} now use theorem \ref{theorem:inv_distribution}} %%\PY{n}{theta\PYZus{}inv\PYZus{}var} \PY{o}{=} \PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{Omega}\PY{p}{,}\PY{n}{Unun}\PY{p}{)}\PY{p}{,}\PY{n}{theta\PYZus{}i\PYZus{}deriv}\PY{o}{.}\PY{n}{transpose}\PY{p}{(}\PY{p}{)}\PY{p}{)} %%\PY{n}{simplify}\PY{p}{(}\PY{n}{theta\PYZus{}inv\PYZus{}var}\PY{p}{)} %%\end{Verbatim} %%\texttt{\color{outcolor}Out[{\color{outcolor}5}]:} %%\begin{equation*} %%\left[\begin{matrix}\frac{2 \mu^{2}}{\sigma^{4}} \left(\mu^{2} + 2 \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{2 \mu^{2}}{\sigma^{4}}\\- \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{1}{\sigma^{4}} \left(2 \mu^{2} + \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}}\\\frac{2 \mu^{2}}{\sigma^{4}} & - \frac{2 \mu}{\sigma^{4}} & \frac{2}{\sigma^{4}}\end{matrix}\right] %%\end{equation*} %%This matches the computation in \eqnref{gauss_confirm_itheta}. Check %%the conjecture: %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}6}]:} \PY{c}{\PYZsh{} now conjecture \ref{conjecture:theta_asym_var_gaussian}} %%\PY{n}{conjec} \PY{o}{=} \PY{n}{quad\PYZus{}form}\PY{p}{(}\PY{n}{Theta\PYZus{}Theta}\PY{p}{,}\PY{n}{Dupp}\PY{p}{)} %%\PY{n}{convar} \PY{o}{=} \PY{l+m+mi}{2} \PY{o}{*} \PY{p}{(}\PY{n}{conjec}\PY{o}{.}\PY{n}{inv}\PY{p}{(}\PY{p}{)} \PY{o}{\PYZhy{}} \PY{n}{Matrix}\PY{p}{(}\PY{l+m+mi}{3}\PY{p}{,}\PY{l+m+mi}{3}\PY{p}{,}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{,}\PY{l+m+mi}{0}\PY{p}{]}\PY{p}{)}\PY{p}{)} %%\PY{n}{simplify}\PY{p}{(}\PY{n}{convar}\PY{p}{)} %%\end{Verbatim} %%\texttt{\color{outcolor}Out[{\color{outcolor}6}]:} %%\begin{equation*} %%\left[\begin{matrix}\frac{2 \mu^{2}}{\sigma^{4}} \left(\mu^{2} + 2 \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{2 \mu^{2}}{\sigma^{4}}\\- \frac{2 \mu}{\sigma^{4}} \left(\mu^{2} + \sigma^{2}\right) & \frac{1}{\sigma^{4}} \left(2 \mu^{2} + \sigma^{2}\right) & - \frac{2 \mu}{\sigma^{4}}\\\frac{2 \mu^{2}}{\sigma^{4}} & - \frac{2 \mu}{\sigma^{4}} & \frac{2}{\sigma^{4}}\end{matrix}\right] %%\end{equation*} %%\begin{Verbatim}[commandchars=\\\{\}] %%{\color{incolor}In [{\color{incolor}7}]:} \PY{c}{\PYZsh{} are they the same?} %%\PY{n}{simplify}\PY{p}{(}\PY{n}{theta\PYZus{}inv\PYZus{}var} \PY{o}{\PYZhy{}} \PY{n}{convar}\PY{p}{)} %%\end{Verbatim} %%\texttt{\color{outcolor}Out[{\color{outcolor}7}]:} %%\begin{equation*} %%\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] %%\end{equation*} %%%UNFOLD %\end{example}%UNFOLD %%UNFOLD \end{document} %for vim modeline: (do not edit) % vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=%%s:syn=rnoweb:ft=rnoweb:nu