{ From user Lonnie, Model Multivariate_distrib at Tue, Nov 11, 2008 1:59 PM~~
}
Softwareversion 4.2.0
Linklibrary Multivariate_distrib
Title: Multivariate Distributions
Description: A library of multivariate distributions.~
~
In a multivariate distribution, each sample is a vector. This vector~~
is identified by an index, identified by the I parameter of the func~~
tions in this library. A Mid value from a distribution function will~~
therefore be indexed by I, whlie a Sample from a distribution functi~~
on is indexed by both I and Run. These distribution functions can al~~
so be used from within the Random function to generate a single monte~~
-carlo sample, which will be indexed by I.~
~
This library also contains functions for generating correlated distri~~
butions. Correlate_with, for example, allows you to generate a univa~~
rite distribution with an arbitrary marginal distribution that has a ~~
specified rank correlation with an arbitrary reference distribution. ~~
Several functions may be used for generating serial correlations, w~~
here each distribution along an index is correlated with the previous~~
point along that index.
Author: Lonnie Chrisman, Ph.D.~
Lumina Decision Systems~
~
With contributions by:~
John Bowers, US FDA.~
Max Henrion, Lumina Decision Systems
Date: Fri, Aug 01, 2003 7:12 PM
Saveauthor: Lonnie
Savedate: Tue, Nov 11, 2008 1:59 PM
Defaultsize: 48,24
Nodesize: 56,24
Nodeinfo: 1,1,1,1,1,1,0,0,0,0
Diagstate: 1,42,10,649,1009,17
Windstate: 2,401,199,483,316
Fontstyle: Arial, 15
Fileinfo: 0,Linklibrary Multivariate_distrib,2,2,0,0,W:\Analytica\Exec~~
Debug\Libraries\Multivariate Distributions.ana
Function Wishart( cv : Number[I,J,Run] ; n :positive ; I,J : Index ; ~
singleSampleMethod : optional hidden scalar)
Title: Wishart(cv,n,I,J)
Description: Suppose you sample N samples from a Gaussian(0,cv,I,J) di~~
stribution, X[I,R]. (R is the index that indexes each sample, R:=1..~~
N). The Wishart distribution describes the distribution of sum( X * ~~
X[I=J], R ). This matrix is dimensioned by I and J and is called the~~
scatter matrix. ~
~
A sample drawn from the Wishart is therefore a sample scatter matrix.~~
If you divide that sample by (N-1), you have a sampled covariance m~~
atrix. ~
~
If you compute a sample covariance matrix from data, and then want to~~
use this in your model, if you just use it directly, you'll be ignor~~
ing sampling error. That may be insignificant of N is large. Otherw~~
ise, you may want to use:~
Wishart( SampleCV, N, I, J) / (N-1)~
instead of just SampleCV in your model. The extended variance will ~~
account for the uncertainty from the finite sample size that was used~~
to obtain your sample CV.~
~
If you can express a prior probability on covariances in the form of ~~
an InvertedWishart distribution, then the posterior distribution, aft~~
er having computed the sample covariance matrix (assumed to be drawn,~~
by nature, from a Wishart), is also an InvertedWishart.
Definition: var T := if i0~
Each sample of a Dirichlet distribution produces a random vector whos~~
e elements sum to 1. It is commonly used to represent second order p~~
robability information.~
~
The Dirichlet distribution has a density given by ~
k * Product( X^(alpha-1), I)~
where k is a normalization factor equal to~
GammaFn( sum(alpha,I )) / Sum(GammaFn(alpha),I)~
~
The parameters, alpha, can be interpreted as observation counts. The~~
mean is given by the relative values of alpha (normalized to 1), but~~
the variance narrows as the alphas get larger, just as your confiden~~
ce in a distribution would narrow as you get more samples.~
~
The Dirichlet lends itself to easy Bayesian updating. If you have a ~~
prior of alpha0, and you observe N
Definition: var a:=Gamma(alpha,singleSampleMethod:singleSampleMethod);~~
~
a/sum(a,I)
Nodelocation: 272,120,1
Nodesize: 58,16
Windstate: 2,26,18,624,485
Paramnames: alpha,I,Over
Function Binormal(MeanVec :numeric[I,Run]; Sdeviations : positive[I,Ru~~
n]; I:IndexType; correlationCoef : numeric[Run];~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar)
Title: BiNormal (m, s, i, c )
Description: A 2-D Normal (or Bi-variate Gaussian) distribution with t~~
he indicated individual standard deviations (>0) and the indicated co~~
rrelation coefficient. The index, I, must have exactly 2 elements, S~~
deviations must be indexed by I.
Definition: if size(I)<>2 then ~
Error("Index to BiNormal must have 2 elements")~
else begin~
var s := product(Sdeviations,I) * correlationCoef;~
Index J:=CopyIndex(I);~
Gaussian( meanVec, If I<>J Then s else Sdeviations^2, I,J,~
singleSampleMethod: singleSampleMethod )~
end
Nodelocation: 288,72,1
Nodesize: 78,16
Windstate: 2,2,24,525,540
Paramnames: MeanVec,Sdeviations,I,correlationCoef,Over
Function Multinomial(N:NonNegative ; theta:NonNegative ; I : IndexTyp~~
e;~
Over : ... optional atomic ;~
singleSampleMethod : hidden optional scalar )
Title: Multinomial (n, theta, i )
Description: Returns the Multinomial Distribution.~
~
The multinomial distribution is a generalization of the Binomial dist~~
ribution to N possible outcomes. For example, if you were to roll a ~~
fair die N times, an outcome would be the number of times each of the~~
six numbers appears. Theta would be the probability of each outcome~~
, where sum(theta,I)=1, and index I is the list of possible outcome. ~~
If theta doesn't sum to 1, it is normalized.~
~
Each sample is a vector indexed by I indicating the number of times t~~
he corresponding outcome (die number) occurred during that sample poi~~
nt. Each sample will have the property that sum( result, I ) = N.
Definition: var z := n;~
var k := size(I);~
~
var j:=cumulate(1,I) in I do begin~
Index I2 := j..k;~
var theta2 := Slice(theta,I,I2); /* unnormalized sub-process */~
var p := theta2/sum(theta2, I2);~
p := if IsNan(p) then 0 else p;~
var xj := Binomial(z,p[I2=j],~
singleSampleMethod:singleSampleMethod);~~
~
z := z - xj;~
xj~
end~
Nodelocation: 117,120,1
Nodesize: 85,16
Windstate: 2,75,167,476,522
Paramnames: N,theta,I,Over
Function Correlate_dists(dists : Context[I,RunIndex] ; rankcorrs : num~~
eric array[I,J] ; ~
I,J : IndexType;~
RunIndex : optional Index = Run )
Title: Correlate Dists (d, rc, i, j )
Description: Reorders the samples in dists so as to match the desired~~
rank correlations between distributions as closely as possible. Ran~~
kCorrs must be positive definite, and the diagonal should contain all~~
ones.~
~
The result will be distributions having the same margins as the origi~~
nal input, but with rank correlations close to those of the rankcorrs~~
matrix.
Definition: if not IsSampleEvalMode and Handle(RunIndex)=Handle(Run) T~~
hen~
dists {Mid mode}~
Else begin~
var u := if Handle(RunIndex)=Handle(run) ~
Then Sample(Gaussian(0,rankcorrs,I,J))~
Else Random(Gaussian(0,rankcorrs,I,J),Over:RunIndex);~
var dsort := sortIndex(dists,RunIndex);~
var urank := UniqueRank(u,RunIndex);~
dists[RunIndex=dsort[RunIndex=urank]]~
end
Nodelocation: 136,392,1
Nodesize: 100,16
Windstate: 2,301,193,557,477
Paramnames: dists,rankcorrs,I,J,RunIndex
Function Correlate_with( S, ref : Context[RunIndex] ; rc : scalar ; ~
RunIndex : optional Index = Run )
Title: Correlate With (s, ref, rc )
Description: Reorders the samples of S so that the result is correlate~~
d with the reference sample with a rank correlation close to rankcorr~~
. ~
~
Example: To generate a logNormal distribution that is highly correlat~~
ed with Ch1, use, e.g.,: Correlate_With( LogNormal(2,3), Ch1, 0.8 )~
~
Note: This achieves a given unweighted rank correlation. If you have~~
a non-default SampleWeighting of points, the weighted rank correlato~~
n may differ.
Definition: if IsSampleEvalMode or Handle(runIndex)<>Handle(Run) Then ~~
begin~
Index q := 1..2;~
var u := If Handle(RunIndex)=Handle(Run) ~
Then binormal( 0, 1, q, rc )~
Else Random(binormal(0,1,q,rc),Over:RunIndex);~
var rrank := UniqueRank(ref,RunIndex);~
var u1sort := sortIndex(u[q=1],RunIndex);~
var u2rank := UniqueRank(u[q=2],RunIndex);~
var ssort := sortIndex(S,RunIndex);~
S[RunIndex=ssort[RunIndex=u2rank[RunIndex=~
u1sort[Ru~~
nIndex=rrank]]]]~
end ~
else {mid mode}~
S
Nodelocation: 128,312,1
Nodesize: 96,16
Windstate: 2,205,170,545,485
Paramnames: S,ref,rc,RunIndex
Function Uniformspherical(I : IndexType ; R : optional Numeric[I,Run] ~~
;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Uniform Spherical (i, r )
Description: Generates points uniformly on a sphere (or circle or hype~~
rsphere).~
Each sample generated is indexed by I -- so if I has 3 elements, the ~~
points will lie on a sphere.~
~
The mid value is a bit strange here since there isn't really a median~~
that lies on the sphere. Obviously the center of the sphere is the ~~
middle value, but that isn't in the allowable range. So, an arbitrar~~
y point on the sphere is used.
Definition: if IsNotSpecified(R) then R:=1;~
var u := Normal(0,1,over:I,~
singleSampleMethod:singleSampleMethod); ~
var d := sqrt( sum(u^2,I) );~
ifall d=0 and @I then R/sqrt(size(I)) else r*u/d
Nodelocation: 328,168,1
Nodesize: 86,16
Windstate: 2,151,227,476,424
Paramnames: I,R,Over
Function Multiuniform(corr : Numeric[I,J,Run] ; I,J : IndexType ; lb,u~~
b : optional Numeric[I,J,Run] ;~
Over : ... optional atomic ;~
singleSampleMethod : hidden optional scalar )
Title: MultiUniform ( c, i, j, lb, ub )
Description: The multi-variate uniform distribution.~
Generates vector samples (indexed by I) such that each component has ~~
a uniform marginal distribution, and such that each component have th~~
e pair-wise correlations given by corr. Indexes I and J must have th~~
e same number of elements, corr needs to be symmetric and must obey a~~
certain semidefinite condition (namely that the transformed matrix [~~
2*sin(30*cov) ] is positive semidefinite. In most cases, this rough~~
ly the same as corr being, or not being, positive semidefinite). Lb ~~
and ub can be used to specify upper and lower bounds, either for all ~~
components, or individually if these bounds are indexed by I. If lb ~~
& ub are omitted, each component will have marginal Uniform(0,1).~
~
The correlation specified in corr is true sample correlation - not ra~~
nk correlation. ~
~
The transformation here is based on:~
* Falk, M. (1999), "A simple approach to the generation of uniformly ~~
distributed random variables with prescribed correlations," Comm. in ~~
Stats - Simulation and Computation 28: 785-791.
Definition: if IsNotSpecified(lb) then lb:=0;~
if IsNotSpecified(ub) then ub := 1;~
var R := if I=J then 1 else 2*sin(30*corr);~
var g := Gaussian(0,R,I,J,~
singleSampleMethod:singleSampleMethod);~~
~
Cumnormal( g ) * (ub-lb) + lb
Nodelocation: 132,168,1
Nodesize: 100,16
Windstate: 2,67,106,608,611
Paramnames: corr,I,J,lb,ub,Over
Module Depricated_multi_var
Title: Depricated multi-variate stuff
Description: Functions found in this module are here for legacy reason~~
s. They existed in older versions of the Multivariate library, but h~~
ave been become obsolete for whatever reason.
Author: Lonnie
Date: Mon, Apr 30, 2007 3:49 PM
Defaultsize: 48,24
Nodelocation: 80,944,1
Nodesize: 56,32
Function Samplecovariance(X ; I : Index ; J : optional Index ; R : Ind~~
ex)
Title: Sample Covariance
Description: This function is obsolete. In Analytica 4.0, the builtin~~
function Variance can be used to compute a covariance matrix. The e~~
quivalent of this function would be: Variance( X, R, CoVarDim:I, CoV~~
arDim2:J ).~
~
Returns a covariance matrix based on the sampled data, X, indexed by ~~
I and R. (I is the dimensionality of X, R corresponds to the samples~~
). The result will be indexed by I and J -- supply J to be the same ~~
length as I.~
~
Note that the mean is simply Average(X,R), and doen't warrant a separ~~
ate function.
Definition: var I2 := if IsNotSpecified(J) ~
Then (Index K/((identifier of I)&"2") := I do VarTerm(K~~
)) ~
Else VarTerm(J);~
var Z:=X-Average(X,R);~
var Zt := Z[@I=@I2];~
Sum(Z*Zt,R)/(size(R)-1)
Nodelocation: 80,48,1
Nodesize: 48,24
Windstate: 2,222,299,476,297
Paramnames: X,I,J,R
Function Samplecorrelation(X : array[I,R] ; I,J,R : IndexType)
Title: sample correlation
Description: This function is obsolete. A covariance matrix can be co~~
mputed in Analytica 4.0+ using the built-in function Correlation. Th~~
e equivalent of this function is Correlation(X,X[@I=@J],R).~
~
Returns a correlation matrix based on data in X, where each data poin~~
t is a vector indexed by I, and the entries in the correlation matrix~~
are the pair-wise correlations of the columns of data. A second ind~~
ex, J, of size identical to I, is required in order to index the 2-di~~
mensional result.
Definition: var z:=x-average(x,R);~
var zt := slice(z,I,cumulate(1,J));~
sum(z*zt,R) / sqrt(sum(z^2,R) * sum(zt^2,R))~
Nodelocation: 208,48,1
Nodesize: 48,24
Windstate: 2,70,24,523,377
Paramnames: X,I,J,R
Close Depricated_multi_var
Text Multvar_te1
Description: Parametric Multivariate Distributions
Nodelocation: 160,40,-1
Nodesize: 136,12
Text Multvar_te2
Description: Creating an array of mutually correlated distributions:
Nodelocation: 232,368,-1
Nodesize: 200,16
Text Multvar_te3
Description: Creating a single univariate distribution correlated wit~~
h another existing dist:
Nodelocation: 296,280,-1
Nodesize: 268,12
Function Normal_correl(m, s, r, y: Numeric ;~
over : optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Normal_correl(m, s, r, y)
Description: Generates a normal distribution with mean m, standard dev~~
iation s, and correlation r with normally distributed value y. In a ~~
deterministic context, it will return m.~
~
If y is not normally distributed, the result will also not be normal,~~
and the correlation will be approximate. It generalizes appropriatel~~
y if any of the parameters are arrays:The result array will have the ~~
union of the indexes of the parameters.
Definition: IF r<-1 OR r>1 THEN Error('Correlation parameter r in func~~
tion Normal_correl(m, s, r, y) is outside the expected range [-1, 1].~~
');~
IFOnly IsSampleEvalMode ~
THEN m + s * (Sqrt(1-r^2) ~
* Normal(Sameindexes( 0, m ), Sameindexes( 1,~~
s ),~
singleSampleMethod:singleSampl~~
eMethod ) ~
+ r * (y - Mean(y))/Sdeviation(y))~
ELSE m
Nodelocation: 352,312,1
Nodesize: 108,16
Windstate: 2,102,90,503,416
Paramnames: m,s,r,y,over
Module Multivariate_interna
Title: Multivariate Internal Functions
Author: Lonnie
Date: Tue, May 01, 2007 9:29 PM
Defaultsize: 48,24
Nodelocation: 200,944,1
Nodesize: 52,32
Diagstate: 1,605,145,550,300,17
Function Sameindexes(x, y)
Title: SameIndexes(x,y)
Description: Returns an array with the same indexes as y, and value x ~~
in each cell.
Definition: IF y=y THEN x ELSE x
Nodelocation: 120,64,1
Nodesize: 80,20
Paramnames: x,y
Function Uniquerank(X : Array[I]; I : Index)
Title: UniqueRank
Description: Returns the Rank of X along I, but such that the rank ass~~
igned is unique for every element. Thus, when there are ties, instea~~
d of getting the same rank, as would happen with the Rank(X,I) functi~~
on, the ranks will be assigned arbitrarily. Consider:~
[ 3, 1, 3, 2, 3, 2, 1 ]~
Ranks become:~
[5,1,6,3,7,4,2 ]
Definition: index Pos := @I;~
var s := SortIndex(X[@I=@Pos],Pos);~
var result := 1;~
for n:=Pos do ( result[@I=s[@Pos=n]] := n );~
result
Nodelocation: 272,64,1
Nodesize: 52,20
Windstate: 2,477,347,537,379
Paramnames: X,I
Close Multivariate_interna
Function Multinormal(m, s: Numeric; cm: ArrayType[i, j,Run]; i , j: In~~
dexType ;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Multinormal(m,s,c,i,j)
Description: A multi-variate normal (or Gaussian) distribution with me~~
an m, standard deviation s, and correlation matrix cm. m and s may ~~
be scalar or indexed by i. cm must be symmetric, positive-definite, a~~
nd indexed by i & j, which must be the same length.~
~
Multinormal uses a correlation matrix. Compare with Gaussian, which ~~
also defines a multi-variate normal but which uses a covariance matri~~
x.
Definition: Gaussian(m,cm*s*s[@i=@j],i,j,over,singleSampleMethod)
Nodelocation: 472,72,1
Nodesize: 84,16
Windstate: 2,391,248,512,343
Paramnames: m,s,cm,i,j,Over
Text Multvar_te4
Description: Reshaped distributions:
Nodelocation: 136,448,-1
Nodesize: 100,16
Function Dist_reshape(x : Numeric[R] ; newdist : all Numeric[R] ; ~
R : optional Index = Run )
Title: Dist_reshape(x, newdist)
Description: Reshapes the probability distribution of uncertain quanti~~
ty x so that it has the same marginal probability distribution (i.e, ~~
same set of sample values) as newdist, but retains the same ranks as ~~
x. Thus:~
Rank(Sample(x), Run) ~
= Rank(Sample(Reshape_dist(x, y)), Run)~
In a Mid context, it simply returns the mid value of newdist, with an~~
y indexes of x.~
~
The result retains any rank correlations that x may have with other p~~
redecessor variables. So, the rank-order correlation between a third~~
variable z and x will be the same as the rank-order correlation betw~~
een z and a reshaped version of x, i.e.~
RankCorrel(x, z) = RankCorrel(Reshape_Dist(x, y), z)~
~
The operation may optionally be applied along an index other than Run~~
.
Definition: IFOnly IsSampleEvalMode or Handle(R)<>Handle(Run) THEN BEG~~
IN~
VAR dsort := SortIndex(newdist, Run);~
VAR xranks := Rank(x, Run);~
newdist[Run = dsort[Run=xranks]]~
END~
ELSE newdist * (x=x)
Nodelocation: 152,472,1
Nodesize: 116,16
Windstate: 2,102,90,646,469
Paramnames: x,newdist,R
Text Multvar_te5
Description: Arrays with serial correlation
Nodelocation: 208,532,-1
Nodesize: 168,12
Function Normal_serial_correl(m, s, r: Numeric; i: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Normal_serial_correl(m,s,r,i)
Description: Generates an array over index i of normal distributions w~~
ith mean m, standard deviation s, and correlation r between successiv~~
e values over index i. You can give each distribution a different m~~
ean and/or standard deviation if m and/or s are arrays indexed by i. ~~
If r is indexed by i, r[i=k] specifies the correlation between result~~
[i=k] and result[i=k-1]. (Then the first correlation, slice(r, i, 1)~~
is ignored.)
Definition: Var x := Normal(0, 1,singleSampleMethod:singleSampleMethod~~
);~
(FOR j := i DO ~
x := Normal_correl( 0, 1, r[i = j],x,~
singleSampleMethod:singleSampl~~
eMethod ) ) ~
* s + m
Nodelocation: 160,560,1
Nodesize: 120,16
Windstate: 2,353,325,540,383
Paramnames: m,s,r,i,over
Function Normal_additive_gro(x, m, s, r: Numeric; i: IndexType ;~
over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Normal_additive_gro(x,m,s,r,i)
Description: Adds a normally distributed percent growth g with mean m ~~
and standard deviation s to x for each value of index i. The growth g~~
for each i has serial correlation r with g for i-1.
Definition: x *( 1 + Cumulate(Normal_serial_correl(m, s, r, i,~
singleSampleMethod:singleSampleMethod), i~~
))
Nodelocation: 159,600,1
Nodesize: 119,16
Windstate: 2,102,90,519,306
Paramnames: x,m,s,r,i,over
Function Normal_compound_gro(x, m, s, r: Numeric; t: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Normal_compound_gro(x,m,s,r,t)
Description: An array of values over time index t, starting from with ~~
value x, and with compound growth applied for each time interval, wit~~
h normal uncertainty with mean m and standard deviation s The growth~~
g for each i has correlation r with g for i-1.
Definition: x * Cumproduct(IF t = Slice(t, 1) THEN 1 ELSE Normal_seria~~
l_correl(m, s, r, t, singleSampleMethod:singleSampleMethod ) + 1, t)~~
Nodelocation: 159,640,1
Nodesize: 119,16
Windstate: 2,102,90,529,366
Paramnames: x,m,s,r,t,over
Function Dist_serial_correl(x; r; i: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Dist_serial_correl(x,r,i)
Description: Generates an array y over index i where each y[i] has a m~~
arginal distribution identical to x, and serial rank correlation of ~~
r with y[i-1]. If x is indexed by i, each y[i] has the same margin~~
al distribution as x[i], but with samples reordered to have the speci~~
fied rank correlation r between successive values. If r is indexed b~~
y i, r[i=k] specifies the rank correlation between y[i=k] and y[i=k-1~~
]. Then the first correlation, r[i=1], is ignored.~
~
In Mid context, it returns Mid(x).~
~
Note: The result retains no probabilistic dependence on x.
Definition: Dist_reshape(Normal_serial_correl( 0, 1, r, i, singleSampl~~
eMethod:singleSampleMethod ), x)
Nodelocation: 408,560,1
Nodesize: 120,16
Windstate: 2,302,78,477,447
Paramnames: x,r,i,over
Function Dist_additive_growth(x, g, r: Numeric; i: IndexType;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Dist_additive_growth(x,g,r,i)
Description: Generates an array of values over index i, with the first~~
equal to x, and successive values adding an uncertain growth with pr~~
obability distribution g, and serial correlation r between growth[i =~~
k] and growth[i=k-1]. x, g, and r each may be indexed by i if you w~~
ant them to vary over i.
Definition: x + Cumulate(Dist_serial_correl( g, r, i, singleSampleMeth~~
od : singleSampleMethod), i)
Nodelocation: 407,600,1
Nodesize: 119,16
Windstate: 2,102,90,506,300
Paramnames: x,g,r,i,over
Function Dist_compound_growth(x, g, r; i: IndexType ;~
over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Dist_compound_growth(x,g,r,i)
Description: Starts with x and applies a compound growth g for each va~~
lue of index i. The growth g for each i has correlation r with g for ~~
i-1.
Definition: x * Cumproduct(~
IF i = Slice(i, 1) THEN 1 ~
ELSE (Dist_serial_correl( g, r, i, ~
singleSampleMethod:singleSampleMe~~
thod ) + 1)~
, i)
Nodelocation: 407,640,1
Nodesize: 119,16
Windstate: 2,102,90,489,307
Paramnames: x,g,r,i,over
Text Multvar_te6
Description: Distributions on Linear Regression coefficients
Nodelocation: 296,688,-1
Nodesize: 256,12
Function Regressionnoise( Y : Numeric[I,Run] ; B : Numeric[I,K,Run] ; ~~
I,K : Index; C : optional Numeric[K,Run] )
Title: RegressionNoise(Y,B,I,K,C)
Description: When you have data, Y[I] and B[I,K], generated from an un~~
derlying model with unknown coefficients C[k] and S of the form:~
~
Y = Sum( C*B, I) + Normal(0,S)~
~
This function computes an estimate for S. ~
~
When using in conjunction with RegressionDist, it is most efficient t~~
o provide the optional parameter C to both routines, where C is the e~~
xpected value of the regression coefficients, obtained from calling R~~
egression(Y,B,I,K). Doing so avoids an unnecessary call to the built~~
in Regression function.
Definition: if IsNotSpecified(C) Then C := Regression(Y,B,I,K);~
Var resid := Y - Sum(C*B,K);~
sqrt( Sum(resid^2,I) / (size(I)-size(K)) );~
Nodelocation: 384,736,1
Nodesize: 104,20
Windstate: 2,332,211,498,542
Paramnames: Y,B,I,K,C
Function Regressionfitprob( Y : Numeric[I,Run] ; B : Numeric[I,K,Run] ~~
; I,K : Index; C : optional Numeric[K,Run] ; ~
S : optional Numeric[I,Run] )
Title: RegressionFitProb(Y,B,I,K,C)
Description: Once you've obtained regression coefficients C (indexed b~~
y K) by calling the Regression function, this function returns the pr~~
obability that a fit this poor would occur by chance, given the assum~~
ption that the data was generated by a process of the form:~
~
Y = Sum( C*B,K) + Normal(0,S)~
~
If this result is very close to zero, it probably indicates that the ~~
assumption of linearity is bad. If it is very close to one, then it ~~
validates the assumption of linearity.~
~
This is not a distribution function - it does not return a sample whe~~
n evaluated in Sample mode. However, it does complement the multivar~~
iate RegressionDist function also included in this library.~
~
To use, first call the Regression function, then you must either know~~
the measurement knows a priori, or obtain it using the RegressionNoi~~
se function.~
~
Var E_C := Regression(Y,B,I,K);~
Var S := RegressionNoise(Y,B,I,K,C);~
Var PrThisPoor := RegressionFitProb(Y,B,I,K,E_C,S)
Definition: if IsNotSpecified(C) then C:=Regression(Y,B,I,K);~
if IsNotSpecified(S) then S:=RegressionNoise(Y,B,I,K);~
var resid := Y - sum(C*B,K);~
var n := size(I);~
var chi2 := sum( resid^2 / Mean(S)^2, I);~
GammaI( n/2 - 1, chi2/2 )
Nodelocation: 152,800,1
Nodesize: 112,20
Windstate: 2,287,69,586,548
Paramnames: Y,B,I,K,C,S
Close Multivariate_distrib