%\VignetteEngine{knitr::knitr} %\VignetteIndexEntry{Using the PDQutils Package} %\VignetteKeyword{distributions} %\VignettePackage{PDQutils} \documentclass[10pt,a4paper,english]{article} % front matter%FOLDUP \usepackage[hyphens]{url} \usepackage{amsmath} \usepackage{amsfonts} % for therefore \usepackage{amssymb} % for theorems? \usepackage{amsthm} % http://en.wikibooks.org/wiki/LaTeX/Rotations \usepackage[counterclockwise]{rotating} \providecommand{\email}[1]{\href{mailto:#1}{\normalfont\texttt{#1}}} \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} \providecommand{\figref}[1]{Figure\nobreakspace\ref{fig:#1}} % 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{PDQutils} %\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{\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=TRUE,cache.path="cache/PDQutils") #opts_chunk$set(fig.path="figure/",dev=c("pdf","cairo_ps")) opts_chunk$set(fig.path="figure/PDQutils",dev=c("pdf")) opts_chunk$set(fig.width="4in",fig.height="4in",fig.width=8,fig.height=8,dpi=450) # 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="small") opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=50,blank=TRUE)) compile.time <- Sys.time() # from the environment # only recompute if FORCE_RECOMPUTE=True w/out case match. FORCE_RECOMPUTE <- (toupper(Sys.getenv('FORCE_RECOMPUTE',unset='False')) == "TRUE") # compiler flags! # not used yet LONG.FORM <- FALSE mc.resolution <- ifelse(LONG.FORM,1000,200) mc.resolution <- max(mc.resolution,100) library(PDQutils) gen_norm <- rnorm lseq <- function(from,to,length.out) { exp(seq(log(from),log(to),length.out = length.out)) } @ %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{Using the \PDQutils package} \author{Steven E. Pav % \thanks{\email{shabbychef@gmail.com}}} %\date{\today, \currenttime} \maketitle %UNFOLD \providecommand{\brakit}[1]{\left[#1\right]} \providecommand{\wrapit}[1]{\left(#1\right)} \providecommand{\ooint}[1]{\left(#1\right)} \providecommand{\coint}[1]{\left[#1\right)} \providecommand{\ocint}[1]{\left(#1\right]} \providecommand{\ccint}[1]{\left[#1\right]} \providecommand{\funcit}[2]{#1\wrapit{#2}} \providecommand{\wt}[1]{\funcit{w}{#1}} \providecommand{\fpi}[2][n]{\funcit{p_{#1}}{#2}} \providecommand{\fexp}[1]{\operatorname{exp}\wrapit{#1}} \providecommand{\He}[2][i]{\funcit{He_{#1}}{#2}} \providecommand{\glag}[3]{\funcit{L^{\wrapit{#1}}_{#2}}{#3}} \providecommand{\jacb}[3]{\funcit{P^{\wrapit{#1}}_{#2}}{#3}} \providecommand{\legn}[2]{\funcit{P_{#1}}{#2}} \providecommand{\chebo}[2]{\funcit{T_{#1}}{#2}} \providecommand{\chebt}[2]{\funcit{U_{#1}}{#2}} \providecommand{\Gam}[1]{\funcit{\Gamma}{#1}} \providecommand{\Bet}[1]{\funcit{B}{#1}} \providecommand{\dx}[1][x]{\mathrm{d}#1} \providecommand{\kdel}[1]{\delta_{#1}} \providecommand{\knn}[1]{h_{#1}} \providecommand{\dens}[1]{\funcit{f}{#1}} \providecommand{\cdf}[1]{\funcit{F}{#1}} \providecommand{\dnorm}[1]{\funcit{\phi}{#1}} \providecommand{\pnorm}[1]{\funcit{\Phi}{#1}} % cf http://people.math.sfu.ca/~cbm/aands/page_775.htm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{abstract}%FOLDUP Example computations via the \PDQutils package are illustrated. \end{abstract}%UNFOLD The \PDQutils package provides tools for approximating the density, distribution, and quantile functions, and for generation of random variates of distributions whose cumulants and moments can be computed. The PDF and CDF are computed approximately via the Gram Charlier A series, while the quantile is computed via the Cornish Fisher approximation. \cite{1998A&AS..130..193B,Jaschke01} The random generation function uses the quantile function and draws from the uniform distribution. \section{Gram Charlier Expansion}%FOLDUP Given the raw moments of a probability distribution, we can approximate the probability density function, or the cumulative distribution function, via a Gram-Charlier A expansion. This is typically developed as an approximation to the normal distribution using Hermite polynomials, but here we follow a more general derivation, which allows us to approximate distributions which are more like a gamma or beta. Let $\wt{x}$ be some non-negative `weighting function', typically the PDF of a known probability distribution. Let $\fpi{x}$ be polynomials which are orthogonal with respect to this weighting function. That is \begin{equation} \int_{-\infty}^{\infty} \wt{x} \fpi[n]{x} \fpi[m]{x} \dx = \kdel{n,m}\knn{n}, \end{equation} where $\kdel{m,n}$ is the Kronecker delta, equal to one only when $m=n$, otherwise equal to zero. We furthermore suppose that the polynomials $\fpi{x}$ are complete: any reasonably smooth function can be represented as a linear combination of these polynomials. Then we can expand the probability density of some random variable, $\dens{x}$ in terms of this basis. Let \begin{equation} \dens{x} = \sum_{n=0}^{\infty} c_n \fpi[n]{x} \wt{x}. \end{equation} By the orthogonality relationship, we can find the constants $c_n$ by multiplying both sides by $\fpi[m]{x}$ and integrating: \begin{equation} \int_{-\infty}^{\infty} \fpi[m]{x} \dens{x} \dx = \sum_{n=0}^{\infty} c_n \int_{-\infty}^{\infty} \fpi[m]{x} \fpi[n]{x} \wt{x} \dx %= \sum_{n=0}^{\infty} c_n \kdel{n,m}\knn{n} = c_m \knn{m}. \end{equation} Thus $$ c_n = \frac{1}{\knn{n}} \int_{-\infty}^{\infty} \fpi[n]{x} \dens{x} \dx $$ When the coefficients of the polynomial $\fpi[n]{x}$ and the uncentered moments of the probability distribution are known, the constant $c_n$ can easily be computed. Thus the density $\dens{x}$ can be approximated by truncating the infinite sum as \begin{equation} \dens{x} \approx \sum_{n=0}^{m} \fpi[n]{x} \wt{x} \brakit{\frac{1}{\knn{n}} \int_{-\infty}^{\infty} \fpi[n]{z} \dens{z} \dx[z]}. \end{equation} To approximately compute the cumulative distribution function, one can compute the integral of the approximate density. The approximation is \begin{equation} \cdf{x} \approx \sum_{n=0}^{m} \int_{-\infty}^{\infty} \fpi[n]{y} \wt{y}\dx[y] \brakit{\frac{1}{\knn{n}} \int_{-\infty}^{\infty} \fpi[n]{z} \dens{z} \dx[z]}. \end{equation} In summary, to approximate the PDF or CDF of a distribution via the Gram Charlier series, one must know the moments of the distribution, and be able to compute $\wt{x}, \fpi[n]{x}, \knn{n},$ and $\int \fpi[n]{y}\wt{y}\dx[y]$. These are collected in Table \ref{tab:orthopoly} for a few different families of probability distributions. \cite[22.2]{abramowitz_stegun} The traditional Gram Charlier `A' series corresponds to case where $\wt{x}$ is the PDF of the standard normal distribution and $\fpi[n]{x}$ is the (probabilist's) Hermite polynomial. Also of interest are the cases where $\wt{x}$ is PDF of the gamma distribution (including Chi-squares), in which case $\fpi[n]{x}$ are the generalized Laguerre polynomials; the case where $\wt{x}$ is the PDF of the (shifted) Beta distribution, and $\fpi[n]{x}$ are the Jacobi polynomials. As special cases of the Beta distribution, one also has the Arcsine distribution (with Chebyshev polynomials of the first kind), the Wigner distribution (Chebyshev of the second kind), and the uniform distribution (Legendre polynomials). \cite{abramowitz_stegun} \nocite{samuelson1943,santos2007,leon2011} %Suppose $f(x)$ is the probability density of some random %variable, and let $F(x)$ be the cumulative distribution function. %Let $He_j(x)$ be the $j$th probabilist's Hermite %polynomial. These polynomials form an orthogonal basis, with respect to the %function $w(x) = e^{-x^2/2} = \sqrt{2\pi}\phi(x)$, of the Hilbert space of %functions which are square integrable with $w$-weighting. \cite[22.2.15]{abramowitz_stegun} %The orthogonality relationship is %$$ %\int_{-\infty}^{\infty} He_i(x) He_j(x) w(x) \mathrm{d}x = \sqrt{2\pi} j! %\delta_{ij}, %$$ %where $\delta_{ij}$ is the Kronecker delta. %Expanding the density $f(x)$ in terms of these polynomials in the %usual way (abusing orthogonality) one has %$$ %f(x) = \sum_{0\le j} \frac{He_j(x)}{j!} \phi(x) \int_{-\infty}^{\infty} f(z) He_j(z) \mathrm{d}z. %$$ %The cumulative distribution function is 'simply' the integral of this %expansion. Abusing certain facts regarding the PDF and CDF of the normal %distribution and the probabilist's Hermite polynomials, the CDF has %the representation %$$ %F(x) = \Phi(x) - \sum_{1\le j} \frac{He_{j-1}(x)}{j!} \phi(x) \int_{-\infty}^{\infty} f(z) He_j(z) \mathrm{d}z. %$$ %These series contain coefficients defined by the probability distribution %under consideration. They take the form %$$ %c_j = \frac{1}{j!}\int_{-\infty}^{\infty} f(z) He_j(z) \mathrm{d}z. %$$ %Using linearity of the integral, these coefficients are easily computed in %terms of the coefficients of the Hermite polynomials and the raw, uncentered %moments of the probability distribution under consideration. Note that it may be the %case that the computation of these coefficients suffers from bad numerical %cancellation for some distributions, and that an alternative formulation %may be more numerically robust. %\begin{sidewaystable}[h] %\centering %\begin{tabular}[b]{|*{6}{c|}} %\hline %class & support & $\wt{x}$ & $\fpi[n]{x}$ & $\knn{n}$ & $\int \fpi[n]{y}\wt{y}\dx[y]$ \\ %\hline %Normal & $\ooint{-\infty,\infty}$ & %$\frac{1}{\sqrt{2\pi}}\fexp{-\frac{x^2}{2}}$ & (Hermite) $\He[n]{x}$ & $n!$ & %$-\frac{1}{\sqrt{2\pi}}\fexp{-\frac{y^2}{2}} \He[n-1]{y}$ \\ %\hline %Gamma & $\coint{0,\infty}$ & $\frac{1}{\Gam{a+1}} x^{a} \fexp{-x}$ & %(generalized Laguerre) $\glag{a}{n}{x}$ & %$\frac{\Gam{n+a+1}}{\Gam{a+1}n!}$ & %$\wrapit{\frac{a+1}{n}} \wrapit{\frac{y^{a+1}\fexp{-y}}{\Gam{a+2}}} \glag{a+1}{n-1}{y}$.\\ %\hline %Beta & $\ccint{-1,1}$ & $\frac{1}{\Bet{a,b}2^{a+b-1}} \wrapit{1-x}^a \wrapit{1+x}^b$ & %(Jacobi) $\jacb{a,b}{n}{x}$ & %$\frac{1}{2n+a+b+1} \frac{\Gam{a+b+2}}{\Gam{a+1}\Gam{b+1}}\frac{\Gam{n+a+1}\Gam{n+b+1}}{n! \Gam{n+a+b+1}}$ & %$\wrapit{\frac{-2}{n}}\wrapit{\frac{\wrapit{a+1}\wrapit{b+1}}{\wrapit{a+b+2}\wrapit{a+b+3}}} %\frac{1}{\Bet{a+2,b+2}2^{a+b+3}}\wrapit{1-y}^{a+1}\wrapit{1+y}^{b+1}\jacb{a+1,b+1}{n-1}{y}$.\\ %\hline %\end{tabular} %\label{table:orthopoly} %\caption{blah blah blah.} %\end{sidewaystable} % % % \begin{table}[htb] \centering \begin{tabular}[b]{|r|*{2}{c|}} \hline class & support & $\wt{x}$ \\ \hline Normal/Hermite & $\ooint{-\infty,\infty}$ & $\dnorm{x} = \frac{1}{\sqrt{2\pi}}\fexp{-\frac{x^2}{2}}$ \\ \hline Gamma/generalized Laguerre & $\coint{0,\infty}$ & $\funcit{g_{a+1}}{x} = \frac{1}{\Gam{a+1}} x^{a} \fexp{-x}$ \\ \hline Beta/Jacobi& $\ooint{-1,1}$ & $\funcit{f_{a+1,b+1}}{x} = \frac{\wrapit{1-x}^a \wrapit{1+x}^b}{\Bet{a+1,b+1}2^{a+b+1}}$ \\ \hline %Arcsine/Chebyshev I & $\ooint{-1,1}$ & $\frac{1}{\pi \sqrt{1-x^2}}$ \\ %\hline %Wigner/Chebyshev II & $\ccint{-1,1}$ & $\frac{2}{\pi} \sqrt{1-x^2}$ \\ %\hline %Uniform/Legendre & $\ccint{-1,1}$ & $\frac{1}{2}$ \\ %\hline \end{tabular} \begin{tabular}[b]{|r|*{2}{c|}} \hline class & $\fpi[n]{x}$ & $\knn{n}$ \\ \hline Normal & $\He[n]{x}$ & $n!$ \\ \hline Gamma & $\glag{a}{n}{x}$ & $\frac{\Gam{n+a+1}}{\Gam{a+1}n!}$ \\ \hline Beta & $\jacb{a,b}{n}{x}$ & $\frac{1}{2n+a+b+1} \frac{1}{\Bet{a+1,b+1}}\frac{\Gam{n+a+1}\Gam{n+b+1}}{n! \Gam{n+a+b+1}}$ \\ \hline %Arcsine & $\chebo{n}{x}$ & $\frac{1 + \kdel{0,n}}{2}$ \\ %\hline %Wigner & $\chebt{n}{x}$ & $1$ \\ %\hline %Uniform & $\legn{n}{x}$ & $\frac{1}{2n+1}$ \\ %\hline \end{tabular} \begin{tabular}[b]{|r|*{1}{c|}} \hline class & $\int \fpi[n]{y}\wt{y}\dx[y]$ \\ \hline Normal & $-\dnorm{y} \He[n-1]{y}$ \\ \hline Gamma & %$\wrapit{\frac{a+1}{n}} \wrapit{\frac{y^{a+1}\fexp{-y}}{\Gam{a+2}}} \glag{a+1}{n-1}{y}$.\\ $\wrapit{\frac{a+1}{n}} \funcit{g_{a+2}}{y} \glag{a+1}{n-1}{y}$.\\ \hline Beta & $\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{a+2,b+2}}{\Bet{a+1,b+1}}} \funcit{f_{a+2,b+2}}{y}\jacb{a+1,b+1}{n-1}{y}$.\\ %$\wrapit{\frac{-2}{n}}\wrapit{\frac{\wrapit{a+1}\wrapit{b+1}}{\wrapit{a+b+2}\wrapit{a+b+3}}} %\frac{1}{\Bet{a+2,b+2}2^{a+b+3}}\wrapit{1-y}^{a+1}\wrapit{1+y}^{b+1}\jacb{a+1,b+1}{n-1}{y}$.\\ \hline %Arcsine & %$\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{3/2,3/2}}{\Bet{1/2,1/2}}} %\funcit{f_{3/2,3/2}}{y}\jacb{1/2,1/2}{n-1}{y}$.\\ %\hline %Wigner & %$\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{5/2,5/2}}{\Bet{3/2,3/2}}} %\funcit{f_{5/2,5/2}}{y}\jacb{3/2,3/2}{n-1}{y}$.\\ %\hline %Uniform & %$\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{2,2}}{\Bet{1,1}}} %\funcit{f_{2,2}}{y}\jacb{1,1}{n-1}{y}$.\\ %\hline \end{tabular} \label{tab:orthopoly} \caption{Different classes of orthogonal polynomials are presented. In each case the weight function, $\wt{x}$ is the PDF of a common distribution, while the orthogonal polynomials come from a well known family. The constant $\knn{n}$ is the normalizing constant. The last table gives the integral of the polynomial times the weighting function, a value which is needed for approximating the CDF. Values are given for: the normal PDF, with probabilist's Hermite polynomials; the Gamma PDF, with generalized Laguerre polynomials; the Beta PDF with Jacobi polynomials. As special cases of the latter, one has the Arcsine, Wigner, and Uniform distributions, with Chebyshev and Legendre polynomials.} \end{table} %\begin{sidewaystable}[h] %\centering %\begin{tabular}[b]{|r|*{5}{c|}} %\hline %class & support & $\wt{x}$ %& $\fpi[n]{x}$ & $\knn{n}$ %& $\int \fpi[n]{y}\wt{y}\dx[y]$ \\ %\hline %Normal/Hermite & $\ooint{-\infty,\infty}$ & %$\dnorm{x} = \frac{1}{\sqrt{2\pi}}\fexp{-\frac{x^2}{2}}$ %& $\He[n]{x}$ & $n!$ %& $-\dnorm{y} \He[n-1]{y}$ \\ %\hline %Gamma/generalized Laguerre & $\coint{0,\infty}$ & $\funcit{g_{a+1}}{x} = \frac{1}{\Gam{a+1}} x^{a} %\fexp{-x}$ %& $\glag{a}{n}{x}$ & $\frac{\Gam{n+a+1}}{\Gam{a+1}n!}$ %& $\wrapit{\frac{a+1}{n}} \funcit{g_{a+2}}{y} \glag{a+1}{n-1}{y}$.\\ %\hline %Beta/Jacobi& $\ooint{-1,1}$ & $\funcit{f_{a+1,b+1}}{x} = \frac{\wrapit{1-x}^a \wrapit{1+x}^b}{\Bet{a+1,b+1}2^{a+b+1}}$ %& $\jacb{a,b}{n}{x}$ & $\frac{1}{2n+a+b+1} \frac{1}{\Bet{a+1,b+1}}\frac{\Gam{n+a+1}\Gam{n+b+1}}{n! \Gam{n+a+b+1}}$ %& $\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{a+2,b+2}}{\Bet{a+1,b+1}}} %\funcit{f_{a+2,b+2}}{y}\jacb{a+1,b+1}{n-1}{y}$.\\ %\hline %Arcsine/Chebyshev I & $\ooint{-1,1}$ & $\frac{1}{\pi \sqrt{1-x^2}}$ %& $\chebo{n}{x}$ & $\frac{1 + \kdel{0,n}}{2}$ %& $\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{3/2,3/2}}{\Bet{1/2,1/2}}} %\funcit{f_{3/2,3/2}}{y}\jacb{1/2,1/2}{n-1}{y}$.\\ %\hline %Wigner/Chebyshev II & $\ccint{-1,1}$ & $\frac{2}{\pi} \sqrt{1-x^2}$ %& $\chebt{n}{x}$ & $1$ %& $\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{5/2,5/2}}{\Bet{3/2,3/2}}} %\funcit{f_{5/2,5/2}}{y}\jacb{3/2,3/2}{n-1}{y}$.\\ %\hline %Uniform/Legendre & $\ccint{-1,1}$ & $\frac{1}{2}$ %& $\legn{n}{x}$ & $\frac{1}{2n+1}$ %& $\wrapit{\frac{-2}{n}}\wrapit{\frac{\Bet{2,2}}{\Bet{1,1}}} %\funcit{f_{2,2}}{y}\jacb{1,1}{n-1}{y}$.\\ %\hline %\end{tabular} %\label{tab:orthopoly_II} %\caption{Different classes of orthogonal polynomials are presented. In each %case the weight function, $\wt{x}$ is the PDF of a common distribution, while %the orthogonal polynomials come from a well known family. The constant %$\knn{n}$ is the normalizing constant. The last table gives the integral of %the polynomial times the weighting function, a value which is needed for %approximating the CDF. Values are given for: the normal PDF, with probabilist's %Hermite polynomials; the Gamma PDF, with generalized Laguerre polynomials; the %Beta PDF with Jacobi polynomials. As special cases of the latter, one has the %Arcsine, Wigner, and Uniform distributions, with Chebyshev and Legendre %polynomials.} %\end{sidewaystable} %UNFOLD \section{Edgeworth Expansion}%FOLDUP Another approximation of the probability density and cumulative distribution functions is the Edgeworth Expansions. These are expressed in terms of the cumulants of the distribution, and also include the Hermite polynomials. However, the derivation of the Edgeworth expansion is rather more complicated than of the Gram Charlier expansion. \cite{1998A&AS..130..193B} The Edgeworth series for a zero-mean unit distribution is $$ f(x) = \frac{1}{\sigma}\phi\left(\frac{x}{\sigma}\right)\left[1 + \sum_{1\le s} \sigma^s \sum_{\{k_m\}} He_{s+2r}\left(x/\sigma\right) \prod_{1 \le m \le s} \frac{1}{k_m!}\left(\frac{S_{m+2}}{(m+2)!}\right)^{k_m}\right], $$ where the second sum is over partitions $\{k_m\}$ such that $k_1 + 2k_2 + \ldots + sk_s = s$, where $r=k_1 + k_2 + \ldots + k_s$, and where $S_n = \frac{\kappa_n}{\sigma^{2n-2}}$ is a semi-normalized cumulant. %UNFOLD \section{Cornish Fisher Approximation}%FOLDUP The Cornish Fisher approximation is the Legendre inversion of the Edgeworth expansion of a distribution, but ordered in a way that is convenient when used on the mean of a number of independent draws of a random variable. Suppose $x_1, x_2, \ldots, x_n$ are $n$ independent draws from some probability distribution. Letting $$ X = \frac{1}{\sqrt{n}} \sum_{1 \le i \le n} x_i, $$ the Central Limit Theorem assures us that, assuming finite variance, $$ X \rightarrow \mathcal{N}(\sqrt{n}\mu, \sigma), $$ with convergence in $n$ The Cornish Fisher approximation gives a more detailed picture of the quantiles of $X$, one that is arranged in decreasing powers of $\sqrt{n}$. The quantile function is the function $q(p)$ such that $P\left(X \le q(p)\right) = q(p)$. The Cornish Fisher expansion is $$ q(p) = \sqrt{n}\mu + \sigma \left(z + \sum_{3 \le j} c_j f_j(z)\right), $$ where $z = \Phi^{-1}(p)$ is the normal $p$-quantile, and $c_j$ involves standardized cumulants of the distribution of $x_i$ of order up to $j$. Moreover, the $c_j$ include decreasing powers of $\sqrt{n}$, giving some justification for truncation. When $n=1$, however, the ordering is somewhat arbitrary. %UNFOLD \section{An Example: Sum of Nakagamis}%FOLDUP The Gram Charlier and Cornish Fisher approximations are most convenient when the random variable can be decomposed as the sum of a small number of independent random variables whose cumulants can be computed. For example, suppose $Y = \sum_{1 \le i \le k} \sqrt{X_i / \nu_i}$ where the $X_i$ are independent central chi-square random variables with degrees of freedom $\nu_1,\nu_2,...,\nu_k$. I will call this a `snak' distribution, since each summand follows a Nakagami distribution. We can easily write code that generates variates from this distribution given a vector of the degrees of freedom: <<'rsnak',eval=TRUE,echo=TRUE>>= rsnak <- function(n,dfs) { samples <- Reduce('+',lapply(dfs,function(k) { sqrt(rchisq(n,df=k)/k) })) } @ Let's take one hundred thousand draws from this distribution. A q-q plot of this sample against normality is shown in \figref{testit}. The normal model is fairly decent, although possibly unacceptable in the tails. Using a Cornish Fisher approximation, we can do better. <<'testit',eval=TRUE,echo=TRUE,fig.width=4.5,fig.height=4,fig.scap=NA,fig.cap=paste0("A q-q plot of ",n.samp," draws from a sum of Nakagamis distribution is shown against normality.")>>= n.samp <- 1e5 dfs <- c(8,15,4000,10000) set.seed(18181) # now draw from the distribution rvs <- rsnak(n.samp,dfs) data <- data.frame(draws=rvs) library(ggplot2) mu <- mean(rvs) sigma <- sd(rvs) ph <- ggplot(data, aes(sample = draws)) + stat_qq(distribution=function(p) { qnorm(p,mean=mu,sd=sigma) }) + geom_abline(slope=1,intercept=0,colour='red') + theme(text=element_text(size=8)) + labs(title="Q-Q plot (against normality)") print(ph) @ Using the additivity property of cumulants, we can compute the cumulants of $Y$ easily if we have the cumulants of the $X_i$. These in turn can be computed from the raw moments. The $j$th moment of a chi distribution with $\nu$ degrees of freedom has form $$ 2^{j/2} \frac{\Gamma\left((\nu + j)/2\right)}{\Gamma\left(\nu/2\right)}. $$ The following function computes the cumulants of the `snak' distribution: <<'snakcu',eval=TRUE,echo=TRUE>>= # for the moment2cumulant function: library(PDQutils) # compute the first ord.max raw cumulants of the sum of Nakagami variates snak_cumulants <- function(dfs,ord.max=10) { # first compute the raw moments moms <- lapply(dfs,function(nu) { ords <- 1:ord.max moms <- 2 ^ (ords/2.0) * exp(lgamma((nu+ords)/2) - lgamma(nu/2)) # we are dividing the chi by sqrt the d.f. moms <- moms / (nu ^ (ords/2.0)) moms }) # turn moments into cumulants cumuls <- lapply(moms,moment2cumulant) # sum the cumulants tot.cumul <- Reduce('+',cumuls) return(tot.cumul) } @ We can now trivially implement the `dpq' functions trivially using the Gram-Charlier and Cornish-Fisher approximations, via \PDQutils, as follows: <<'dpqsnak',eval=TRUE,echo=TRUE>>= library(PDQutils) dsnak <- function(x,dfs,ord.max=10,...) { raw.moment <- cumulant2moment(snak_cumulants(dfs,ord.max)) retval <- dapx_gca(x,raw.moment,support=c(0,Inf),...) return(retval) } psnak <- function(q,dfs,ord.max=10,...) { raw.moment <- cumulant2moment(snak_cumulants(dfs,ord.max)) retval <- papx_gca(q,raw.moment,support=c(0,Inf),...) return(retval) } qsnak <- function(p,dfs,ord.max=10,...) { raw.cumul <- snak_cumulants(dfs,ord.max) retval <- qapx_cf(p,raw.cumul,support=c(0,Inf),...) return(retval) } @ An alternative version of the PDF and CDF functions using the Edeworth expanion would look as follows: <<'dpqsnak2',eval=TRUE,echo=TRUE>>= dsnak_2 <- function(x,dfs,ord.max=10,...) { raw.cumul <- snak_cumulants(dfs,ord.max) retval <- dapx_edgeworth(x,raw.cumul,support=c(0,Inf),...) return(retval) } psnak_2 <- function(q,dfs,ord.max=10,...) { raw.cumul <- snak_cumulants(dfs,ord.max) retval <- papx_edgeworth(q,raw.cumul,support=c(0,Inf),...) return(retval) } @ Using this approximate quantile function, the q-q plot looks straighter, as shown in \figref{improvedqq}. <<'improvedqq',eval=TRUE,echo=TRUE,fig.width=4.5,fig.height=4,fig.scap=NA,fig.cap=paste0("A q-q plot of ",n.samp," draws from a sum of Nakagamis distribution is shown against quantiles from the `qsnak' function.")>>= data <- data.frame(draws=rvs) library(ggplot2) ph <- ggplot(data, aes(sample = draws)) + stat_qq(distribution=function(p) { qsnak(p,dfs=dfs) }) + geom_abline(slope=1,intercept=0,colour='red') + theme(text=element_text(size=8)) + labs(title="Q-Q against qsnak (C-F appx.)") print(ph) @ Note that the q-q plot uses the approximate quantile function, computed via the Cornish-Fisher expansion. We can test the Gram Charlier expansion by computing the approximate CDF of the random draws and checking that it is nearly uniform, as shown in \figref{psnak_ecdf}. <<'psnak_ecdf',eval=TRUE,echo=TRUE,fig.width=4.5,fig.height=4,fig.scap=NA,fig.cap=paste0("The empirical CDF of the approximate CDF of a sum of Nakagamis distribution on ",n.samp," draws is shown.")>>= apx.p <- psnak(rvs,dfs=dfs) require(ggplot2) ph <- ggplot(data.frame(pv=apx.p),aes(x=pv)) + stat_ecdf(geom='step') print(ph) @ %UNFOLD \section{A warning on convergence} Blinnikov and Moessner note that the the Gram Charlier expansion will actually diverge for some distributions when more terms in the expansion are considered, behaviour which is not seen for the Edgeworth expansion. \cite{1998A&AS..130..193B} Here, we will replicate their example of the chi-square distribution with 5 degrees of freedom. Blinnikov and Moessner actually transform the chi-square to have zero mean and unit variance. They plot the true PDF of this normalized distribution, along with the 2- and 6-term Gram Charlier approximations, as shown in \figref{chisetup}. <<'chisetup',eval=TRUE,echo=TRUE,fig.width=4.5,fig.height=4,fig.scap=NA,fig.cap="The true PDF of a normalized $\\chi^2_5$ distribution is shown, along with the 2- and 6-term Gram Charlier approximations. This replicates Figure 1 of Blinnikov and Moessner. \\cite{1998A&AS..130..193B}">>= # compute moments and cumulants: df <- 5 max.ord <- 20 subords <- 0:(max.ord - 1) raw.cumulants <- df * (2^subords) * factorial(subords) raw.moments <- cumulant2moment(raw.cumulants) # compute the PDF of the rescaled variable: xvals <- seq(-sqrt(df/2) * 0.99,6,length.out=1001) chivals <- sqrt(2*df) * xvals + df pdf.true <- sqrt(2*df) * dchisq(chivals,df=df) pdf.gca2 <- sqrt(2*df) * dapx_gca(chivals,raw.moments=raw.moments[1:2],support=c(0,Inf)) pdf.gca6 <- sqrt(2*df) * dapx_gca(chivals,raw.moments=raw.moments[1:6],support=c(0,Inf)) all.pdf <- data.frame(x=xvals,true=pdf.true,gca2=pdf.gca2,gca6=pdf.gca6) # plot it by reshaping and ggplot'ing require(reshape2) arr.data <- melt(all.pdf,id.vars='x',variable.name='pdf',value.name='density') require(ggplot2) ph <- ggplot(arr.data,aes(x=x,y=density,group=pdf,colour=pdf)) + geom_line() print(ph) @ Compare this with the 2 and 4 term Edgeworth expansions, shown in \figref{chitwo}. <<'chitwo',eval=TRUE,echo=TRUE,fig.width=4.5,fig.height=4,fig.scap=NA,fig.cap="The true PDF of a normalized $\\chi^2_5$ distribution is shown, along with the 2- and 4-term Edgeworth expansions. This replicates Figure 6 of Blinnikov and Moessner. \\cite{1998A&AS..130..193B}">>= # compute the PDF of the rescaled variable: xvals <- seq(-sqrt(df/2) * 0.99,6,length.out=1001) chivals <- sqrt(2*df) * xvals + df pdf.true <- sqrt(2*df) * dchisq(chivals,df=df) pdf.edgeworth2 <- sqrt(2*df) * dapx_edgeworth(chivals,raw.cumulants=raw.cumulants[1:4],support=c(0,Inf)) pdf.edgeworth4 <- sqrt(2*df) * dapx_edgeworth(chivals,raw.cumulants=raw.cumulants[1:6],support=c(0,Inf)) all.pdf <- data.frame(x=xvals,true=pdf.true,edgeworth2=pdf.edgeworth2,edgeworth4=pdf.edgeworth4) # plot it by reshaping and ggplot'ing require(reshape2) arr.data <- melt(all.pdf,id.vars='x',variable.name='pdf',value.name='density') require(ggplot2) ph <- ggplot(arr.data,aes(x=x,y=density,group=pdf,colour=pdf)) + geom_line() print(ph) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % bibliography%FOLDUP \nocite{edgeworthCF,1998A&AS..130..193B,cheah1993,AS269,Jaschke01,CFetc} \bibliographystyle{plainnat} \bibliography{PDQutils} %UNFOLD \end{document} %for vim modeline: (do not edit) % vim:fdm=marker:fmr=FOLDUP,UNFOLD:cms=%%s:syn=rnoweb:ft=rnoweb:nu