Nothing
# A Poisson Subordinated Distribution By Lihn (LIHNPSD)
library(moments)
library(BB) # spg
library(Bolstad2)
library(optimx)
library(Rmpfr) # Multiple Precision Float, erf, erfc are defined here !
library(sn) # skew normal
.N <- function(x,p=128) mpfr(x, precBits = p)
# library(codetools)
# findGlobals(calcqq)
# ------------------------------------------
# LIHNPSD_UnitTest(): for unit testing
prepare <- function(d) UseMethod("prepare", d)
rawmean <- function(d) UseMethod("rawmean", d)
rawmu1 <- function(d) UseMethod("rawmu1", d)
rawmu2 <- function(d) UseMethod("rawmu2", d)
rawmu3 <- function(d) UseMethod("rawmu3", d)
rawmu4 <- function(d) UseMethod("rawmu4", d)
rawsn <- function(d, type, x, k) UseMethod("rawsn", d)
density <- function(d,x) UseMethod("density", d)
rawdensity <- function(d,x) UseMethod("rawdensity", d)
rawdensity_kth <- function(d,x,k) UseMethod("rawdensity_kth", d)
rawdensityslope <- function(d,x) UseMethod("rawdensityslope", d)
rawdensity0 <- function(d) UseMethod("rawdensity0", d)
rawcdf <- function(d,x) UseMethod("rawcdf", d)
rawcdfinv <- function(d,c,xinit) UseMethod("rawcdfinv", d)
tailindex <- function(d,x) UseMethod("tailindex", d)
tailindex_plot <- function(d,xmin,xmax,ymax=0.4) UseMethod("tailindex_plot", d)
psdkernel <- function(d,k) UseMethod("psdkernel", d)
psdmagnitude <- function(d,r,step=0.1) UseMethod("psdmagnitude", d)
poisson_sum <- function(d,fn) UseMethod("poisson_sum", d)
poisson_sum_kth <- function(d,fn,k) UseMethod("poisson_sum_kth", d)
psdunittest <- function(d) UseMethod("psdunittest", d)
generatepdf <- function(d,NS,NT,raw=1) UseMethod("generatepdf", d)
calcqq <- function(d,hq,step=5, debug=0) UseMethod("calcqq", d)
mu1_analytic <- function(d) UseMethod("mu1_analytic", d)
mu2_analytic <- function(d) UseMethod("mu2_analytic", d)
mu3_analytic <- function(d) UseMethod("mu3_analytic", d)
mu4_analytic <- function(d) UseMethod("mu4_analytic", d)
mu_n_core <- function(d,n) UseMethod("mu_n_core", d)
psdvariance <- function(d) UseMethod("psdvariance", d)
psdskewness <- function(d) UseMethod("psdskewness", d)
psdkurtosis <- function(d) UseMethod("psdkurtosis", d)
standardfit <- function(d, r, hist, trace, iter, plotqq, weights, merge_tails) UseMethod("standardfit", d)
# when this is in a package
#formals(sd.default) <- c(formals(sd.default), alist(... = ))
# ------------------------------------------
# simple PSD constructor
SPSD <- function(sigma, alpha, gamma, beta=0, mpfr=0) {
d <- list( sigma= sigma, alpha= alpha, gamma= gamma, beta= beta, mpfr= mpfr )
class(d) <- "LIHNPSD"
rawmean(d)
}
# the input can be numeric or character
# notice that 0.05 in double is not exactly that in mpfr
prepare.LIHNPSD <- function(d) {
# you have to provide sigma, alpha and gamma
stopifnot( ! is.null(d$sigma) )
stopifnot( ! is.null(d$alpha) )
stopifnot( ! is.null(d$gamma) )
if (is.null(d$mpfr)) d$mpfr <- 0
D <- function(x) {
if ( d$mpfr > 0 ) {
if ( class(x) != "mpfr" ) mpfr(x, precBits = floor(d$mpfr))
else x # already mpfr
}
else as.numeric(x)
}
if (is.null(d$location)) d$location <- 0.0
if (is.null(d$beta)) d$beta <- 0.0
if (is.null(d$lambda)) d$lambda <- 1.0
if (is.null(d$epsilon)) d$epsilon <- 1e-10
d$pi <- if ( d$mpfr > 0 ) Const("pi", prec= d$mpfr) else d$pi <- pi
d$sigma <- D(d$sigma)
d$alpha <- D(d$alpha)
d$gamma <- D(d$gamma)
d$location <- D(d$location)
d$beta <- D(d$beta)
d$lambda <- D(d$lambda)
d$epsilon <- D(d$epsilon)
d
}
# ------------------------------------------
# PoissonTerm <- function(lambda,k) dpois(k,lambda) # doesn't work in MPFR
PoissonTerm <- function(lambda,k) lambda^k*exp(-lambda)/factorial(k)
LihnFunctionKth <- function(alpha,lambda,k) {
PoissonTerm(lambda,k)*(k+1)^alpha
}
LihnFunctionSum <- function(alpha,x,epsilon = 1e-10) {
# s: sum, k: poisson index, v: value of k-th term
s <- k <- v <- 0.0 * alpha # multiply by alpha to inherit its type
while ( k < 20 || v > epsilon*s || v == 0 ) {
v <- LihnFunctionKth(alpha,x,k)
s <- s + v
k <- k + 1
}
s
}
LihnFunctionAnalytic <- function(alpha,x) {
if ( x <= 0 ) stopifnot( paste("Unsupported x:",as.numeric(x)) == "STOP" )
if ( alpha == 1 ) (x+1)
else if ( alpha == 2 ) (x^2+3*x+1)
else if ( alpha == 3 ) (x^3+6*x^2+7*x+1)
else if ( alpha == 4 ) (x^4+10*x^3+25*x^2+15*x+1)
else if ( alpha == 0 ) 1
else if ( alpha == -1 && x != 0 ) (1-exp(-x))/x
else stopifnot( paste("Unsupported alpha:",as.numeric(alpha)) == "STOP" )
}
LihnBetaPoly <- function(N,b) {
if ( N == 1 ) b
else if ( N == 2 ) 1
else if ( N == 4 ) 3
else if ( N == 6 ) 15
else {
ppi <- if ( class(b) != "mpfr" ) pi
else Const("pi", prec= getPrec(b))
if ( N == 3 ) (-ppi/2*b^3+3*b)
else if ( N == 5 ) (3*ppi^2/4*b^5-5*ppi*b^3+15*b)
else stopifnot( paste("Unsupported N:",as.numeric(N)) == "STOP" )
}
}
LihnTildeFunction <- function(x,alpha,p,epsilon = 1e-10) {
# s: sum, k: poisson index, v: value of k-th term
s <- k <- v <- 0.0 * alpha # multiply by alpha to inherit its type
while ( k < 20 || v > epsilon*s || v == 0 ) {
sigma <- (k+1)^alpha
v <- (k+1)^p/factorial(k) * exp(-x^2/(2*sigma^2))
s <- s + v
k <- k + 1
}
C <- if ( class(alpha) == "mpfr" )
exp(-1+0*alpha)/sqrt(2* Const("pi", prec= getPrec(alpha)))
else exp(-1)/sqrt(2*pi)
s * C
}
ParetoDiff2 <- function(x,alpha) {
L5 <- LihnTildeFunction(x,alpha,-5*alpha)
L3 <- LihnTildeFunction(x,alpha,-3*alpha)
A <- (x^2*L5- L3)
L1 <- LihnTildeFunction(x,alpha,-alpha)
L3 <- LihnTildeFunction(x,alpha,-3*alpha);
A/L1 - x^2 * L3^2/L1^2
}
LihnFunctionValidate <- function() {
# we will validate Lihn function at alpha =-1,0,1,2,3,4
err <- 1e-6
for (a in c(-1,0,1,2,3,4)) {
for( x in c(1,2,3,4,5) ) {
# print(paste("LihnFunction: validating for a=",
# as.numeric(a),"x=",as.numeric(x)))
stopifnot( abs(LihnFunctionSum(a,x) - LihnFunctionAnalytic(a,x)) < err )
}
}
print("LihnFunction is validated")
stopifnot( abs( LihnTildeFunction(0,0,0) - 1/sqrt(2*pi) ) < err )
stopifnot( abs( LihnTildeFunction(1,0,0) - 0.24197072451914 ) < err )
stopifnot( abs( LihnTildeFunction(1,0.5,-0.5) - 0.21950993766309 ) < err )
stopifnot( abs( ParetoDiff2(0,0) + 1 ) < err )
stopifnot( abs( ParetoDiff2(1,0) + 1 ) < err )
stopifnot( abs( ParetoDiff2(1,0.5) + 0.57402540354604 ) < err )
stopifnot( abs( ParetoDiff2(10,0.62) + 0.018150532857661 ) < err )
print("LihnTildeFunction is validated");
}
# ------------------------------------------
psdkernel.LIHNPSD <- function(d,k) ((k+1)^d$alpha) * ((1+d$gamma)^k)
# return the k magnitude of the return series
psdmagnitude.LIHNPSD <- function(d, r, step=0.1) {
ks <- c()
for ( i in 1:length(r) ) {
j <- 0
repeat {
if ( abs(r[i]) <= psdkernel(d,j) * d$sigma ) {
ks[i] <- j
break
}
j <- j + step
}
}
ks
}
# This is the utility to build moments
poisson_sum_kth.LIHNPSD <- function(d, fn, k) {
PoissonTerm(d$lambda, k) * fn(d,k)
}
poisson_sum.LIHNPSD <- function(d, fn) {
# sm: sum, m: k-th term of sum
sm <- k <- m <- 0.0*d$sigma # multiply by sigma to inherit its type
kmax <- if ( class(d$sigma) == "mpfr" ) 100000 else 160
cond <- TRUE
while ( cond && k <= kmax ) {
m <- poisson_sum_kth(d, fn, k)
sm <- sm + m
cond <- ( k < 20
|| ( abs(m) >= abs(d$epsilon*sm) && sm != 0.0 )
|| ( sm == 0.0 && k < 20 ) )
# print(paste(cond,"k=",as.numeric(k),"sm=",as.numeric(sm),
# "m=",as.numeric(m),"err=",as.numeric(d$epsilon)))
k <- k + 1
}
sm
}
# ------------------------------------------
rawmu1.LIHNPSD <- function(d) {
mu1core <- function(d,k) ( psdkernel(d, k) * d$sigma )
poisson_sum(d, mu1core) * d$beta
}
rawmean.LIHNPSD <- function(d) {
d <- prepare(d)
d$rawmean <- rawmu1(d)
d
}
rawmu2.LIHNPSD <- function(d) {
mu2core <- function(d,k) {
sigma2 <- psdkernel(d, k) * d$sigma
sigma2^2
}
poisson_sum(d, mu2core)
}
rawmu3.LIHNPSD <- function(d) {
mu3core <- function(d,k) {
sigma2 <- psdkernel(d, k) * d$sigma
sigma2^3
}
poisson_sum(d, mu3core) * LihnBetaPoly(3,d$beta)
}
rawmu4.LIHNPSD <- function(d) {
mu4core <- function(d,k) {
sigma2 <- psdkernel(d, k) * d$sigma
sigma2^4
}
poisson_sum(d, mu4core) * LihnBetaPoly(4,d$beta)
}
# ------------------------------------------
# also raw moments, but expressed in Lihn functions
mu_n_core.LIHNPSD <- function(d, n) {
L <- LihnFunctionSum(n*d$alpha, d$lambda*(1+d$gamma)^n, d$epsilon)
d$sigma^n * exp(d$lambda*((1+d$gamma)^n-1.0)) * L # muN core component
}
mu1_analytic.LIHNPSD <- function(d) {
d <- prepare(d)
L <- LihnFunctionSum(d$alpha, d$lambda*(1+d$gamma), d$epsilon)
mu1 <- d$beta * d$sigma * exp(d$lambda*d$gamma) * L
mu1
}
mu2_analytic.LIHNPSD <- function(d) {
d <- prepare(d)
L <- LihnFunctionSum(2*d$alpha, d$lambda*(1+d$gamma)^2, d$epsilon)
mu2 <- d$sigma^2 * exp(d$lambda*d$gamma*(d$gamma+2)) * L
mu2
}
mu3_analytic.LIHNPSD <- function(d) {
d <- prepare(d)
mu3 <- LihnBetaPoly(3,d$beta) * mu_n_core(d,3)
mu3
}
mu4_analytic.LIHNPSD <- function(d) {
d <- prepare(d)
mu4 <- LihnBetaPoly(4,d$beta) * mu_n_core(d,4)
mu4
}
# ------------------------------------------
psdvariance.LIHNPSD <- function(d) {
rawmu2(d) - d$rawmean^2
}
# http://en.wikipedia.org/wiki/Moment_about_the_mean#Relation_to_moments_about_the_origin
psdskewness.LIHNPSD <- function(d) {
mu1 <- d$rawmean
(rawmu3(d) - 3*mu1*rawmu2(d) + 2*mu1^3 ) / psdvariance(d)^(3/2)
}
psdkurtosis.LIHNPSD <- function(d) {
mu1 <- d$rawmean
(rawmu4(d) - 4*mu1*rawmu3(d) + 6*mu1^2*rawmu2(d) - 3*mu1^4 ) / psdvariance(d)^2 - 3.0
}
# ------------------------------------------
# PDF offset by mu1+d$location
density.LIHNPSD <- function(d, x) {
d2 <- d
d2$location <- 0
rawdensity(d2, x - d$rawmean - d$location)
}
# PDF without any offset (not even d$location)
# k >= 0 (int) is used to study the kth term contribution
rawdensity_kth.LIHNPSD <- function(d, x, k) {
gauss <- function(d,k2) {
sigma2 <- psdkernel(d, k2) * d$sigma
mn <- d$beta * sigma2
gs <- exp(-(x-mn)^2/(2*sigma2^2))/sqrt(2*d$pi)/sigma2
gs
}
poisson_sum_kth(d, gauss, k)
}
# rawsn is an internal function
# type 1: for pdf, 2: for 2nd term of dP/dx
rawsn.LIHNPSD <- function(d, type, x, k) {
sigma2 <- psdkernel(d, k) * d$sigma
delta <- d$beta * sqrt(d$pi/2);
if ( abs(delta) >= 1 ) stopifnot( paste("d$beta is out of range:",as.numeric(d$beta)) == "STOP" )
if ( type == 1 ) {
sqrt2 <- sqrt( d$lambda / d$lambda * 2 ) # accomodate mpfr
sn_a <- delta/sqrt(1-delta^2);
gs <- exp(-x^2/(2*sigma2^2))/sqrt(2*d$pi)/sigma2
gcdf <- (1+erf(sn_a * x/sqrt2/sigma2))/2
gs * 2*gcdf
}
else if ( type == 2 ) {
x2 <- x * sqrt( 2/d$pi ) / sqrt( 2/d$pi - d$beta^2 )
gs <- exp(-x2^2/(2*sigma2^2))/sqrt(2*d$pi)/sigma2
gcdf <- (1+erf(0))/2
gs * 2*gcdf
}
else stopifnot( paste("Unsupported type:",as.numeric(type)) == "STOP" )
}
rawdensity.LIHNPSD <- function(d, x) {
sqrt2 <- sqrt( d$lambda / d$lambda * 2 ) # accomodate mpfr
delta <- d$beta * sqrt(d$pi/2);
if ( abs(delta) >= 1 ) stopifnot( paste("d$beta is out of range:",as.numeric(d$beta)) == "STOP" )
sn_a <- delta/sqrt(1-delta^2);
gauss <- function(d,k) rawsn(d, 1, x, k)
poisson_sum(d, gauss)
}
rawdensityslope.LIHNPSD <- function(d, x) {
dpk_dx <- function(d,k) {
sigma2 <- psdkernel(d, k) * d$sigma
t1 <- -rawsn(d, 1, x, k) * ( x / sigma2^2 )
b2 <- d$beta / sigma2 * sqrt( 2/d$pi ) / sqrt( 2/d$pi - d$beta^2 )
t2 <- b2 * rawsn(d, 2, x, k)
t1 + t2
}
poisson_sum(d, dpk_dx)
}
rawcdf.LIHNPSD <- function(d, x) {
# erf <- function(x) 2*pnorm(x*sqrt(2))-1
# we are using erf from rmpfr package
sqrt2 <- sqrt( d$lambda / d$lambda * 2 ) # accomodate mpfr
delta <- d$beta * sqrt(d$pi/2);
if ( abs(delta) >= 1 ) stopifnot( paste("d$beta is out of range:",as.numeric(d$beta)) == "STOP" )
sn_a <- delta/sqrt(1-delta^2);
cdf <- function(d,k) {
sigma2 <- psdkernel(d, k) * d$sigma
gcdf <- (1+erf(x/sqrt2/sigma2))/2
gcdf - 2*T.Owen(x/sigma2, sn_a)
}
poisson_sum(d, cdf)
}
# xinit: if you don't know what to supply, use d$rawmean
rawcdfinv.LIHNPSD <- function(d, c, xinit) {
err <- d$epsilon
x <- xinit
x1 <- d$rawmean # a fake number
repeat {
pdf <- rawdensity(d,x)
cdf <- rawcdf(d,x)
x1 <- x - (cdf-c)/pdf
if ( abs(x1-x) < err ) break
x <- x1
}
x1
}
rawdensity0.LIHNPSD <- function(d) {
gs <- 1/sqrt(2*d$pi)/d$sigma
f <- exp(-d$lambda*d$gamma/(d$gamma+1))
L <- LihnFunctionSum( -d$alpha, d$lambda/(d$gamma+1) )
gs * f * L
}
tailindex.LIHNPSD <- function(d,x) {
cdf <- rawcdf(d,x)
pdf <- rawdensity(d,x)
dff <- rawdensityslope(d,x)
-1-(1-cdf)*dff/pdf^2
}
tailindex_plot.LIHNPSD <- function(d,xmin,xmax,ymax=0.4) {
N <- 1000
x <- xmin+(0:(N-1))*(xmax-xmin)/N
idx <- tailindex(d,x)
ymax <- min(abs(ymax), max(abs(as.numeric(idx))*2))
plot( x, idx, pch = ".", col="red",
ylim=c(-ymax,+ymax),
ylab="Tail Index(x)", xlab="x", main="Tail Index")
abline(h=0)
# text to describe d
}
# -----------------------------------------
generatepdf.LIHNPSD <- function(d, NS, NT, raw=1) {
d <- prepare(d)
# NS: number of sigma, NT: number of tick sample per unit of sigma
XM <- d$sigma*NS # max of x
N <- NS*NT # number of samples in x-axis
x <- seq(-XM, XM, length=N+1)
dx <- (2*XM)/N
pb <- if ( raw==1 ) rawdensity(d,x) else density(d,x)
list( NS = NS, NT = NT, XM = XM,
N = N, x = x, dx = dx, pb = pb
)
}
# hq = list( qhx = h$mids, qhy = h$counts )
calcqq.LIHNPSD <- function(d, hq, step=5, debug=0 ) {
tm0 <- unclass(Sys.time())
hq$qhyc <- cumsum(hq$qhy) / sum(hq$qhy)
qqp <- list ( x = c(), xq = c(), y = c(), yq = c() )
xmin <- min(hq$qhx)
xmax <- max(hq$qhx)
xabs <- max(abs(xmin),abs(xmax))
tx <- d$location # starting point of theoretical x
dx <- ( hq$qhx[2] - hq$qhx[1] )/step
while ( tx >= hq$qhx[1] ) tx <- tx-dx
stopifnot( hq$qhyc[1] != 0 )
tm1 <- unclass(Sys.time()) - tm0
while( rawcdf(d, tx-d$rawmean-d$location) >= hq$qhyc[1]
&& abs(tx) < 5*xabs ) tx <- tx-dx
if ( debug > 0 ) print(paste("tx=",as.numeric(tx)))
tm2 <- unclass(Sys.time()) - tm0
tcdf <- 0
cdfcnt <- 0
for ( i in 1:(length(hq$qhy)-1) ) {
qqp$x[i] <- hq$qhx[i] # data's x
qqp$xq[i] <- qval <- hq$qhyc[i] # data's cdf(x)
# the goal is to find the min(x) of fit whose cdf >= data's cdf
repeat {
tcdf <- rawcdf(d, tx-d$rawmean-d$location)
cdfcnt <- cdfcnt + 1
if ( tcdf >= qval || abs(tx) > 5*xabs ) break
tx <- tx+dx
}
qqp$y[i] <- tx # fit's x
qqp$yq[i] <- tcdf # fit's cdf(x)
if ( debug == 2 ) print(paste("qx=",as.numeric(qqp$x[i]),
"qval=",as.numeric(qval),"tcdf=",as.numeric(tcdf),"tx=",as.numeric(tx)))
}
if ( debug > 0 ) print(paste("qqp=","tcdf=",as.numeric(tcdf),"tx=",as.numeric(tx)))
tm3 <- unclass(Sys.time()) - tm0
qqp$tm1 <- tm1
qqp$tm2 <- tm2
qqp$tm3 <- tm3
print(paste("qqp tm1=",sprintf("%.2f",tm1),
"tm2=",sprintf("%.2f",tm2),
"tm3=",sprintf("%.2f",tm3),
"cdfcnt=",cdfcnt))
qqp
}
# ------------------------------------------
# Unit Test Functions
# ------------------------------------------
psdunittest.LIHNPSD <- function(d) {
d <- prepare(d)
d <- rawmean(d)
err <- 1e-6
errn <- 1e-6 # for the extreme case
dgs <- 7 # round
Q <- function(x) round(as.numeric(x),dgs) # helper
TM <- function() format(Sys.time(), "%X %x") # helper
print(paste("--Unit Test on: sigma=",as.numeric(d$sigma),
"alpha=",Q(d$alpha),
"gamma=",Q(d$gamma),
"beta=",Q(d$beta),
"lambda=",Q(d$lambda),
"mpfr=",Q(d$mpfr),"--"))
# when symmetric, test a few obvious facts
if ( d$beta == 0 ) {
# cdf at x=0 is 1/2
cdf0 <- rawcdf(d,0)
print(paste("cdf_at_0 =",Q(cdf0)))
stopifnot( abs(cdf0 - 0.5) < 1e-4 )
# d pdf/dx at x=0 is 0
stopifnot( rawdensityslope(d,0) == 0 )
}
# test location (only applicable to density, not rawdensity)
d$location <- d$sigma
stopifnot( abs( density(d,d$sigma) - rawdensity(d,d$sigma-d$rawmean-d$location) ) < err )
d$location <- 0
stopifnot( abs( density(d,d$sigma) - rawdensity(d,d$sigma-d$rawmean) ) < err )
NS <- 20 # number of sigma
repeat {
pdf <- rawdensity(d,d$sigma*NS)
if ( pdf <= 1e-20 || NS >= 400 ) break
NS <- NS+20
if ( NS %% 80 == 0 ) print(paste("Pdf test: # sigma=",NS,"pdf=",as.numeric(pdf),"time=",TM() ))
}
NT <- floor(40000/NS) # number of tick sample per unit of sigma
print(paste("Integral: #sigma=",NS,"#ticks=",NT,"pdf=",as.numeric(pdf), "time=",TM() ))
rs <- generatepdf(d, NS, NT)
N <- rs$N
x <- rs$x
pb <- rs$pb
print(paste("PDF generated at ",TM() ))
mu0_num <- sintegral( x, pb, N/2 )$int
mu1_num <- sintegral( x, x * pb, N/2 )$int
mu2_num <- sintegral( x, x^2 * pb, N/2 )$int
mu3_num <- sintegral( x, x^3 * pb, N/2 )$int
mu4_num <- sintegral( x, x^4 * pb, N/2 )$int
# total of probably density should be one
if ( !( abs(1.0 - mu0_num) < err ) ) {
print(paste("mu0_num=",Q(mu0_num)))
stopifnot( abs(1.0 - mu0_num) < err )
}
# validate cdf, seem a bit off that I need to use err=1e-4
L <- length(x)
for( i in floor(L* c(0.45, 0.5-NT/L/2, 0.5, 0.5+NT/L*2, 0.55, 1.0)) ) {
xm <- x[i]
cdf_raw <- rawcdf(d,xm)
cdf_num <- sintegral( x[1:i], pb[1:i], i/2 )$int
print(paste("cdf[",i,"](",Q(xm),")_raw=",
Q(cdf_raw),"_num=",Q(cdf_num)
))
stopifnot( abs(cdf_raw - cdf_num) < 1e-4 )
# pdf slope testing
slope_num <- (rawdensity(d,x[i])-rawdensity(d,x[i-1]))/(x[i]-x[i-1])
slope_raw <- rawdensityslope(d,x[i])
print(paste("slope[",i,"](",Q(xm),")_raw=",
Q(slope_raw),"_num=",Q(slope_num)
))
stopifnot( abs(slope_raw - slope_num) < 5e-3 ) # very loose
}
pb0_num <- rawdensity(d,0)
pb0_raw <- rawdensity0(d)
print(paste("pb0_raw=",Q(pb0_raw),Q(pb0_num)))
stopifnot( abs(pb0_raw - pb0_num)/pb0_raw < err )
mu2_raw <- rawmu2(d)
mu3_raw <- rawmu3(d)
mu4_raw <- rawmu4(d)
var_raw <- psdvariance(d)
sk_raw <- psdskewness(d)
kt_raw <- psdkurtosis(d)
var_num <- mu2_num - mu1_num^2
sk_num <- (mu3_num - 3*mu1_num*mu2_num + 2*mu1_num^3 ) / var_num^(3/2)
kt_num <- (rawmu4(d) - 4*d$rawmean*rawmu3(d) + 6*d$rawmean^2*rawmu2(d) - 3*d$rawmean^4 ) / var_num^2 - 3.0
print(paste("var_raw=",Q(var_raw),Q(var_num),
"err=", as.numeric(abs(var_raw-var_num)/var_raw)
))
print(paste("sk_raw=",Q(sk_raw),Q(sk_num),
"err=", as.numeric(abs(sk_raw-sk_num)/max(1,as.numeric(abs(sk_raw))))
))
print(paste("kt_raw=",Q(kt_raw),Q(kt_num),
"err=", as.numeric(abs(kt_raw-kt_num)/max(1,as.numeric(abs(kt_raw))))
))
mu1_a <- mu1_analytic(d)
mu2_a <- mu2_analytic(d)
mu3_a <- mu3_analytic(d)
mu4_a <- mu4_analytic(d)
print(paste(
"rawmean=",Q(d$rawmean)
,"mu1_a=",Q(mu1_a)
,"mu1_n=",Q(mu1_num)
))
if ( kt_raw > 10 ) errn <- 1e-4
stopifnot( abs(d$rawmean - mu1_a) < err )
stopifnot( abs(d$rawmean - mu1_num) < errn )
print(paste(
"mu2_raw=",Q(mu2_raw)
,"mu2_a=",Q(mu2_a)
,"mu2_n=",Q(mu2_num)
))
if ( kt_raw > 10 ) errn <- 1e-3
stopifnot( abs(mu2_raw - mu2_a)/mu2_raw < err )
stopifnot( abs(mu2_raw - mu2_num)/mu2_raw < errn )
stopifnot( abs(var_raw - var_num)/var_raw < errn )
print(paste(
"mu3_raw=",Q(mu3_raw)
,"mu3_a=",Q(mu3_a)
,"mu3_n=",Q(mu3_num)
))
if ( kt_raw > 10 ) errn <- 1e-2
dn <- max(abs(mu3_raw),mu2_raw)
stopifnot( abs(mu3_raw - mu3_a)/dn < err )
stopifnot( abs(mu3_raw - mu3_num)/dn < errn )
stopifnot( abs(sk_raw - sk_num) < errn )
print(paste(
"mu4_raw=",Q(mu4_raw)
,"mu4_a=",Q(mu4_a)
,"mu4_n=",Q(mu4_num)
))
if ( kt_raw > 10 ) errn <- 1/20
stopifnot( abs(mu4_raw - mu4_a)/mu4_raw < err )
stopifnot( abs(mu4_raw - mu4_num)/mu4_raw < errn )
stopifnot( abs(kt_raw - kt_num)/max(0.1,as.numeric(kt_raw)) < errn )
# CDF^-1 test
CDFCUT <- 0.01 # 1 percent
xmax <- 0
repeat {
if ( rawcdf(d,xmax-1) < 2*CDFCUT ) break
xmax <- xmax-1
}
xmin <- xmax-1
repeat {
if ( rawcdf(d,xmin-1) < 0.3*CDFCUT ) break
xmin <- xmin-1
}
xi <- seq(xmin,xmax,by=0.001)
ci <- rawcdf(d,xi)
x_num <- xi[max(which(ci<CDFCUT))] # solved by brute force
x_ntn <- rawcdfinv(d, CDFCUT, d$rawmean) # solved by Newton's method
print(paste(
"CDF=",CDFCUT
,"x_num=",Q(x_num)
,"x_ntn=",Q(x_ntn)
))
stopifnot( abs(x_num - x_ntn)/abs(x_ntn) < 1e-3 )
print(paste("----","----",TM(),"----","----"))
}
LIHNPSD_UnitTest <- function(mpfr=0) {
tm <- format(Sys.time(), "%X %x")
print(paste("----","----",tm,"----","----"))
LihnFunctionValidate()
# In mpfr, fractional number needs to be quoted to be precise
# simple case
dist <- list( sigma= 1.0, alpha= 0.0, gamma= 0, beta= "-0.1", mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
dist <- list( sigma= 1.0, alpha= 0.5, gamma= "0.2", beta= 0, mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
dist <- list( sigma= 1.0, alpha= 0.5, gamma= 0, beta= "-0.1", mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
# symmetric case
dist <- list( sigma= 1.0, alpha= 1.0, gamma= 0, beta= 0, mpfr= mpfr)
class(dist) <- "LIHNPSD"
psdunittest(dist)
dist <- list( sigma= 1.0, alpha= 1.0, gamma= 0, beta= -0.25, mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
dist <- list( sigma= 1.0, alpha= 1.0, gamma= "0.3", beta= 0.0, mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
# Approximately DJI
# We need to study why this is not converging?
# var err ~ 8e-6, can't get much better no matter NT, NS, epsilon
dist <- list( sigma= 1.0, alpha= 1.0, gamma= "0.3", beta= "-0.2", mpfr= mpfr )
class(dist) <- "LIHNPSD"
psdunittest(dist)
# Maxima: rawmu1(0.002,2,0.3,-0.1,1,300),numer; -0.00177911390838
print("--LIHNPSD_UnitTest is completed successfully--")
}
# ------------------------------------------
# Fitting Functions
# ------------------------------------------
LIHNPSD_standardfit_test <- function(d, r, hist, plotqq= 1,
weights=list(), merge_tails=c(0,0) ) {
data_stats <- c( mean(r), sqrt(var(r)), skewness(r), kurtosis(r)-3.0 )
psd <- c( d$location/data_stats[2], d$sigma/data_stats[2], d$alpha, d$gamma, d$beta )
LIHNPSD_standardfit_fn(psd, data_stats, hist, weights=weights, plotqq= plotqq, merge_tails= merge_tails)
}
LIHNPSD_standardfit_fn <- function( psd, data_stats, hist, plotqq=1,
weights=list(), merge_tails=c(0,0), debug=0 ) {
tm0 <- unclass(Sys.time())
d <- list( location= psd[1]*data_stats[2],
sigma= psd[2]*data_stats[2], alpha= psd[3], gamma= psd[4], beta= psd[5] )
class(d) <- "LIHNPSD"
d <- rawmean(d)
if ( abs(d$beta) >= sqrt(2/pi) ) {
delta <- if ( debug == 1 ) list( delta= abs(d$beta)*1000 ) else abs(d$beta)*1000
return(delta)
}
if (is.null(weights$m2)) weights$m2 <- 1
if (is.null(weights$m3)) weights$m3 <- 1
if (is.null(weights$m4)) weights$m4 <- 1
if (is.null(weights$mmx)) weights$mmx <- 1
if (is.null(weights$pdf_df)) weights$pdf_df <- 1
if (is.null(weights$qq_df)) weights$qq_df <- 0.1
if (is.null(weights$a_bump)) weights$a_bump <- 1
st <- c( 0, sqrt(psdvariance(d)), psdskewness(d), psdkurtosis(d) )
tm1 <- unclass(Sys.time()) - tm0
NS <- 4
xs <- seq( d$rawmean-NS*d$sigma, d$rawmean+NS*d$sigma, by=2*NS*d$sigma/200 )
ds <-rawdensity(d,xs)
max_density <- max(hist$density)
mmx <- abs( max(ds)-max_density ) / max_density
tm2 <- unclass(Sys.time()) - tm0
NS <- 4
pdf_df <- 0.0
xpdf <- c()
tpdf <- c()
j <- 1
for ( i in 1:length(hist$mids) ) {
x <- hist$mids[i]
if ( abs(x-d$rawmean-d$location) <= NS*d$sigma ) {
xpdf[j] <- x
tpdf[j] <- rawdensity(d, x-d$rawmean-d$location)
df <- (log(hist$density[i]) - log(tpdf[j]))^2
pdf_df <- pdf_df + if ( abs(df) <= 1000.0 ) df else 0
j <- j + 1
}
}
pdf_df <- sqrt(abs(pdf_df))
pdf_df <- if ( pdf_df == Inf ) 10000.0 else pdf_df
tm3 <- unclass(Sys.time()) - tm0
hq <- MergeTailHistogram( list(qhx=hist$mids, qhy=hist$counts), merge_tails)
qqp <- calcqq(d, hq, debug=debug)
qq_df <- sqrt(sum(abs(qqp$y-qqp$x)^2))/max(hist$mids)
tm4 <- unclass(Sys.time()) - tm0
alpha_bump <- if ( d$alpha > 1.5 ) abs(d$alpha-1.5)/10 else 0
m2 <- abs( st[2]-data_stats[2] ) / data_stats[2] # var pct
m3 <- abs( st[3]-data_stats[3] )
m4 <- abs( st[4]-data_stats[4] ) / data_stats[4] # kurt pct
delta <- 0
delta <- delta + (m2 * weights$m2)^2
delta <- delta + (m3 * weights$m3)^2
delta <- delta + (m4 * weights$m4)^2
delta <- delta + (alpha_bump * weights$a_bump)^2
delta <- delta + (mmx * weights$mmx)^2
delta <- delta + (qq_df * weights$qq_df)^2
delta <- delta + (pdf_df * weights$pdf_df)^2
delta <- sqrt(delta)
if ( plotqq == 1 ) {
par(mfcol=c(2,2)); # 1 plot
par(oma = c(1, 1, 4, 1))
par(mar=c(4,4,2,2))
xmin <- min(xs+d$rawmean+d$location)
xmax <- max(xs+d$rawmean+d$location)
plot(hist$mids, hist$density, type="o", pch=22, lty=2, col="red",
main="PDF", xlab="r", ylab="PDF",
xlim=c(xmin,xmax) )
lines(xpdf, tpdf, pch = ".", col="black")
lines(xs+d$rawmean+d$location, ds, pch = ".", col="green")
par(mar=c(4,4,2,2))
plot(hist$mids, log(hist$density), type="o", pch=22, lty=2, col="red",
main="Log PDF", xlab="log(r)", ylab="Log PDF")
lines(xpdf, log(tpdf), pch = ".", col="black")
xmin <- min(hist$mids)
xpos <- xmin*0.5
ypos <- max(log(hist$density))
sc <- 0.8
text(xmin*0.5,ypos,labels=c("data"),cex=sc)
text(xmin*0.5,ypos-1,labels=paste("mean", sprintf("%.6f",d$rawmean)),cex=sc)
text(xmin*0.5,ypos-2,labels=paste("std", sprintf("%.6f",sqrt(psdvariance(d)))),cex=sc)
text(xmin*0.5,ypos-3,labels=paste("skew", sprintf("%.6f",psdskewness(d))),cex=sc)
text(xmin*0.5,ypos-4,labels=paste("kurt", sprintf("%.6f",psdkurtosis(d))),cex=sc)
par(mar=c(4,4,2,2))
plot(qqp$x, qqp$y, pch = "o", col="red",
main=paste("QQ Plot",Sys.time()),
ylim=c( min(qqp$y), max(qqp$y)*1.5 ),
xlab="Observed Quantile", ylab="Theoretical Quantile")
abline(h=0,col="yellow")
abline(v=0,col="yellow")
lines(qqp$x, abs(qqp$y-qqp$x), pch = ".", col="green")
abline(0,1,col="blue")
tm9 <- unclass(Sys.time()) - tm0
xmin <- min(qqp$x)
xpos <- xmin*0.5
ypos <- max(qqp$y)*1.5
ydf <- (max(qqp$y)-min(qqp$y))/10
sc <- 0.8
text(xpos,ypos-1*ydf,labels=c("timing"),cex=sc)
text(xpos,ypos-2*ydf,labels=paste("tm1", sprintf("%.2f",tm1)),cex=sc)
text(xpos,ypos-3*ydf,labels=paste("tm2", sprintf("%.2f",tm2)),cex=sc)
text(xpos,ypos-4*ydf,labels=paste("tm3", sprintf("%.2f",tm3)),cex=sc)
text(xpos,ypos-5*ydf,labels=paste("tm4", sprintf("%.2f",tm4)),cex=sc)
text(xpos,ypos-6*ydf,labels=paste("tm9", sprintf("%.2f",tm9)),cex=sc)
par(mar=c(4,4,2,2))
plot(qqp$x, qqp$y, pch = "o", col="green",
main=paste("QQ Data",Sys.time()),
ylim=c( min(qqp$y), max(qqp$y)*1.5 ),
xlab="Observed Quantile", ylab="Theoretical Quantile")
text(xpos,ypos-1*ydf,labels=c("fit"),cex=sc)
text(xpos,ypos-2*ydf,labels=paste("delta", sprintf("%.8f",delta)),cex=sc)
text(xpos,ypos-4*ydf,labels=paste("alpha", sprintf("%.8f",d$alpha)),cex=sc)
text(xpos,ypos-5*ydf,labels=paste("gamma", sprintf("%.8f",d$gamma)),cex=sc)
text(xpos,ypos-6*ydf,labels=paste("beta", sprintf("%.8f",d$beta)),cex=sc)
text(xpos,ypos-7*ydf,labels=paste("sigma", sprintf("%.6f %.8f",d$sigma,psd[2])),cex=sc)
text(xpos,ypos-8*ydf,labels=paste("loc", sprintf("%.6f %.8f",d$location,psd[1])),cex=sc)
xmax <- max(qqp$x)
xpos <- xmax*0.4
text(xpos,ypos-1*ydf,labels=c("delta parts"),cex=sc)
text(xpos,ypos-2*ydf,labels=paste("max_pdf", sprintf("%.6f",max_density)),cex=sc)
text(xpos,ypos-3*ydf,labels=paste("mmx", sprintf("%.6f %.2f",mmx, weights$mmx)),cex=sc)
text(xpos,ypos-4*ydf,labels=paste("m2", sprintf("%.6f %.2f",m2, weights$m2)),cex=sc)
text(xpos,ypos-5*ydf,labels=paste("m3", sprintf("%.6f %.2f",m3, weights$m3)),cex=sc)
text(xpos,ypos-6*ydf,labels=paste("m4", sprintf("%.6f %.2f",m4, weights$m4)),cex=sc)
text(xpos,ypos-7*ydf,labels=paste("a_bump", sprintf("%.6f %.2f",alpha_bump, weights$a_bump)),cex=sc)
text(xpos,ypos-8*ydf,labels=paste("pdf_df", sprintf("%.6f %.2f",pdf_df, weights$pdf_df)),cex=sc)
text(xpos,ypos-9*ydf,labels=paste("qq_df", sprintf("%.8f %.2f",qq_df, weights$qq_df)),cex=sc)
}
if ( debug == 1 ) list( delta= delta,
mmx= mmx, m2=m2, m3=m3, m4=m4, alpha_bump=alpha_bump,
max_density=c(max(ds),max_density), qq_df= qq_df )
else delta
}
# optimx: http://cran.r-project.org/web/packages/optimx/optimx.pdf
# spg: http://cran.r-project.org/web/packages/BB/BB.pdf
standardfit.LIHNPSD <- function(d, r, hist, trace= 0, iter=5000,
plotqq= 1, weights=list(), merge_tails=c(0,0) ) {
print(paste("standardfit:","iter=",iter,"plotqq=",plotqq,
"weights=(",paste(weights),")"))
data_stats <- c( mean(r), sqrt(var(r)), skewness(r), kurtosis(r)-3.0 )
psd <- c( d$location/data_stats[2], d$sigma/data_stats[2], d$alpha, d$gamma, d$beta )
psdout <- optimx( psd, LIHNPSD_standardfit_fn, method=c("spg"),
itnmax=iter, control=list(trace=trace,eps=1e-5),
data_stats= data_stats, hist= hist, plotqq= plotqq,
weights=weights, merge_tails=merge_tails )
par <- psdout$par$par
d2 <- list( location=par[1]*data_stats[2],
sigma=par[2]*sqrt(var(r)), alpha=par[3], gamma=par[4], beta=par[5] )
class(d2) <- "LIHNPSD"
d2 <- rawmean(d2)
list( dist= d2, psdout= psdout)
}
# ------------------------------------------
# Misc Helper Functions
# ------------------------------------------
# converts the time series to log returns
# to calculates the log-return distribution
TimeSeriesLogReturn <- function( pr, days ) {
n <- length(pr)
s <- seq(1, n, by=days)
pr2 <- pr[s]
r <- c(0.0)
for (i in 2:length(pr2)) {
if ( pr[i-1] > 0.0 && pr[i] ) {
r[i-1] <- log(pr[i]/pr[i-1])
}
}
r
}
ConvertDistToSample <- function( x, hx, N, ngc = 200 ) {
NR <- N / sum(hx)
hx_r <- c()
for( i in 1:length(hx) ) {
xi <- x[i]
xs <- seq(xi,xi,length=hx[i]*NR)
hx_r <- c(hx_r, xs)
if ( i %% ngc ) gc()
}
gc()
hx_r
}
# q = list( qhx = h$mids, qhy = h$counts )
MergeTailHistogramOneSide <- function( q, allowed_merge ) {
qhx <- q$qhx
qhy <- q$qhy
owed <- 0
while ( allowed_merge > 0 ) {
if ( qhy[1] > 0 ) {
qhy <- qhy[-1]
qhx <- qhx[-1]
allowed_merge <- allowed_merge - 1
owed <- owed + 1
while ( qhy[1] == 0 ) {
qhy <- qhy[-1]
qhx <- qhx[-1]
}
}
else break
}
qhy[1] <- qhy[1] + owed
list( qhx= qhx, qhy= qhy )
}
MergeTailHistogram <- function ( q, merge_tails ) {
allowed_negative_merge <- merge_tails[1]
allowed_positive_merge <- merge_tails[2]
q2 <- MergeTailHistogramOneSide ( q, allowed_negative_merge )
q2 <- list( qhx= rev(q2$qhx), qhy= rev(q2$qhy) )
q3 <- MergeTailHistogramOneSide ( q2, allowed_positive_merge )
list( qhx= rev(q3$qhx), qhy= rev(q3$qhy) )
}
# ------------------------------------------------------------------
# Standard plot facility
# ------------------------------------------------------------------
# h: histogram
# tx: theoretical x, tpdf/tcdf: theoretical PSD
LIHNPSD_plotpdf <- function (dist, h, tx, tpdf, xlab="log(r)", main="PSD PDF" ) {
par(mar=c(4,4,2,1))
s <- sqrt(psdvariance(dist))
i <- 3
repeat {
if ( rawcdf(dist, dist$rawmean+i*s) >= 0.995 || i > 20 ) break
i <- i + 1
}
xmin <- -s*i+dist$rawmean
ymax <- max(max(tpdf),max(h$density))
plot( h$mids, h$density, pch = "o", col="red",
xlim=c(-s*i, s*i)+dist$rawmean,
ylim=c(0,ymax),
ylab="PDF", xlab=xlab, main=main)
lines(tx, tpdf, pch = ".", col="blue")
sc=0.8
yseg <- 0.20*max(h$density)
segments( dist$location+dist$rawmean-s, yseg,
dist$location+dist$rawmean+s, yseg, lwd=3 )
text( dist$location, yseg*0.5, "std dev", cex=sc )
# legends
legend(xmin, ymax, c("data","fit"), cex=sc,
col=c("red","blue"), pch=c("o",NA), lty=c(NA,1));
}
LIHNPSD_plotlogpdf <- function (dist, h, tx, tpdf, xlab="log(r)", main="PSD Log PDF" ) {
par(mar=c(4,4,2,1))
fnt <- log(h$density) [ which(is.finite(log(h$density))) ]
xmin <- min(h$mids)
xmax <- max(h$mids)
plot( h$mids, log(h$density), pch = "o", col="red",
ylim=c( min(max(fnt)-8,min(fnt)), max(fnt)+0.5 ),
xlim=c( xmin, xmax ),
ylab="log(PDF)", xlab=xlab, main=main)
lines(tx, log(tpdf), pch = ".", col="blue")
# legends
# parameters
xpos <- 0.7 * if ( abs(xmin) > abs(xmax) ) xmin else xmax
ypos <- max(log(h$density))+0.5
sc <- 0.8
sp <- 0.8
text(xpos,ypos-1*sp,labels=c("psd fit"),cex=sc)
text(xpos,ypos-2*sp,labels=paste("location", sprintf("%.6f",dist$location)),cex=sc)
text(xpos,ypos-3*sp,labels=paste("sigma", sprintf("%.6f",dist$sigma)),cex=sc)
text(xpos,ypos-4*sp,labels=paste("alpha", sprintf("%.4f",dist$alpha)),cex=sc)
text(xpos,ypos-5*sp,labels=paste("gamma", sprintf("%.4f",dist$gamma)),cex=sc)
text(xpos,ypos-6*sp,labels=paste("beta", sprintf("%.4f",dist$beta)),cex=sc)
}
LIHNPSD_plotcdf <- function (dist, h, data_st, tx, tcdf, xlab="log(r)", main="PSD CDF" ) {
par(mar=c(4,4,2,1))
s <- sqrt(psdvariance(dist))
i <- 3
repeat {
if ( rawcdf(dist, dist$rawmean+i*s) >= 0.995 || i > 20 ) break
i <- i + 1
}
plot( h$mids, cumsum(h$counts) / sum(h$counts), lwd=2, type="s", col="red",
xlim=c(-s*i, s*i)+dist$rawmean,
ylab="CDF", xlab=xlab, main=main)
dx <- tx[2]-tx[1]
# lines( h$mids, cumsum(h$counts) / sum(h$counts), col="red")
lines(tx, tcdf, pch = ".", col="blue")
xmin <- -s*i+dist$rawmean
xmax <- s*i+dist$rawmean
sc <- 0.8
text(xmin*0.5,0.90,labels=c("data"),cex=sc)
text(xmin*0.5,0.75,labels=paste("mean", sprintf("%.6f",data_st[1])),cex=sc)
text(xmin*0.5,0.60,labels=paste("std", sprintf("%.4f",data_st[2])),cex=sc)
text(xmin*0.5,0.45,labels=paste("skew", sprintf("%.4f",data_st[3])),cex=sc)
text(xmin*0.5,0.30,labels=paste("kurt", sprintf("%.4f",data_st[4])),cex=sc)
text(xmax*0.5,0.75,labels=c("psd fit"),cex=sc)
text(xmax*0.5,0.60,labels=paste("mean", sprintf("%.6f",dist$rawmean+dist$location)),cex=sc)
text(xmax*0.5,0.45,labels=paste("std", sprintf("%.4f",sqrt(psdvariance(dist)))),cex=sc)
text(xmax*0.5,0.30,labels=paste("skew", sprintf("%.4f",psdskewness(dist))),cex=sc)
text(xmax*0.5,0.15,labels=paste("kurt", sprintf("%.4f",psdkurtosis(dist))),cex=sc)
}
LIHNPSD_plotqq <- function (dist, qqp, merge_tails, main="PSD QQ-Plot" ) {
par(mar=c(4,4,2,1))
dx <- qqp$x[2] - qqp$x[1]
plot(qqp$x+dx/2, qqp$y, pch="o", col="red", main=main,
xlab="Observed Quantile", ylab="Theoretical Quantile")
abline(0,1,col="blue")
abline(h=0,col="yellow")
abline(v=0,col="yellow")
lines(qqp$x, abs(qqp$y-qqp$x), pch = ".", col="green")
# legends
sc <- 0.8
xmin <- min(qqp$x)
ymax <- max(qqp$y)
legend(xmin, ymax, c("qq data","45 degree","error"), cex=sc,
col=c("red","black","green"), pch=c("o",NA,NA), lty=c(NA,1,1));
xmax <- max(qqp$x)
ymin <- min(qqp$y)
text(xmax*0.5, ymin*0.4, "tail dropped:",cex=sc);
text(xmax*0.5, ymin*0.6, paste("left", merge_tails[1]),cex=sc);
text(xmax*0.5, ymin*0.8, paste("right", merge_tails[2]),cex=sc);
}
LIHNPSD_plot_std4gr <- function ( th, dt, EPS=FALSE, file=NA ) {
if ( EPS ) {
postscript(file=file,
paper="special", width=8, height=8, horizontal=FALSE)
}
par(mfcol=c(2,2)); # 4 plots
par(oma = c(1, 1, 4, 1))
LIHNPSD_plotpdf (th$dist, dt$h, th$tx, th$tpdf)
LIHNPSD_plotlogpdf (th$dist, dt$h, th$tx, th$tpdf)
LIHNPSD_plotcdf (th$dist, dt$h, dt$stats, th$tx, th$tcdf)
LIHNPSD_plotqq (th$dist, th$qqp, th$merge_tails )
if ( EPS ) dev.off()
}
LIHNPSD_prepare_data <- function( logr, breaks, merge_tails ) {
dt <- list( logr= logr, N= length(logr), breaks=breaks, merge_tails=merge_tails )
dt$stats <- c( mean(logr), sqrt(var(logr)), skewness(logr), kurtosis(logr)-3.0 )
h <- hist(logr, breaks = breaks, plot = FALSE)
h$cuml <- cumsum(h$counts) / sum(h$counts)
dt$hq <- MergeTailHistogram( list(qhx=h$mids, qhy=h$counts), merge_tails )
dt$h <- h
dt
}
LIHNPSD_theoretical_result <- function ( dist, dt, N=5000 ) {
th <- list( dist=dist, N= N, qqp= calcqq(dist, dt$hq), merge_tails=dt$merge_tails )
tx <- seq(min(dt$logr), max(dt$logr), length=N+1)
th$dx <- tx[2]-tx[1]
th$tpdf <- density(dist,tx)
th$tcdf <- rawcdf(dist, tx-dist$rawmean-dist$location)
th$tx <- tx
#k <- seq(0,10,by=1)
#kn <- psdkernel(dist, k)
#kp <- PoissonTerm( dist$lambda, k)
th
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.