Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.soft-sys.math.mathematica > #3660 > unrolled thread

Re: MultinormalDistribution Question

Back to article view | Back to comp.soft-sys.math.mathematica

This discussion starts older than the indexed window; earlier articles aren't shown. The article labeled Started by below is the oldest one visible, not the original post.

Contents

  Re: MultinormalDistribution Question Andrzej Kozlowski <akoz@mimuw.edu.pl> - 2011-07-11 11:01 +0000
    Re: MultinormalDistribution Question Steve <s123@epix.net> - 2011-07-17 10:03 +0000
      Re: MultinormalDistribution Question Ray Koopman <koopman@sfu.ca> - 2011-07-19 10:54 +0000
      Re: MultinormalDistribution Question Ray Koopman <koopman@sfu.ca> - 2011-07-20 10:36 +0000

#3660 — Re: MultinormalDistribution Question

First, your covariance matrix is not symmetric definite. I think you probably meant:

CapSigma = {{sigma1^2, rho*sigma1*sigma2}, {rho*sigma1*sigma2,
   sigma2^2}}

Assuming that, you can get the conditional expectation from the definition:

Integrate[y PDF[dist, {6.3, y}], {y, -Infinity, Infinity}]/
 PDF[MarginalDistribution[dist, 1], 6.3]

5.6

Alternatively, you can use Mathematica 8 built in NExpectation function:

NExpectation[y \[Conditioned] 6.299 <= x <= 6.301, {x, y} \[Distributed] dist]
5.6


Once you have the conditional expectation, you can compute the conditional variance, e.g.

 Chop[Integrate[(y - 5.6)^2*PDF[dist, {6.3, y}], {y, -Infinity, Infinity}]/PDF[MarginalDistribution[dist, 1], 6.3]]
 0.0256


Andrzej Kozlowski

On 10 Jul 2011, at 11:01, Steve wrote:

> Hello,
>
> Can someone help me with this ?
>
> I have 2 normal distributions; dist1 describes x and dist2 describes
> y. Each are fully defined and are correlated to one another by the
> correlation coefficient. How can I detemine the mean and standard
> deviation of the expected normal distribution that is associated with
> a given x value from dist1 ?
>
> An example:
> mean1 = 5.8
> sigma1 =0 .2
>
> mean2 = 5.3
> sigma2 = 0.2
>
> Correlation Coefficient, rho = 0.6
>
> Given an x value of 6.3 (from dist1) what is the corresponding mean
> and standard deviation of y ?
>
> I can view the combined density function from the following:
>
> Mu = {mean1, mean2}
> CapSigma = {{sigma1^2, rho*sigma1*sigma2} , {rho,  rho*sigma1*sigma2}
> dist = MultinormalDistribution[Mu,CapSigma]
> pdf = PDF[dist,{x,y}]
> plot1 = Plot3D[pdf, {x,4,7},{y,4,7}, PlotRange->All]
>
> but can't see how to determine the mean and the standard deviation of
> y for a given value of x, like 6.3
>
> Any help would be appreciated.
>
> Thanks,
>
> --Steve
>
>
>
>

[toc] | [next] | [standalone]

#3776

On Jul 11, 7:01 am, Andrzej Kozlowski <a...@mimuw.edu.pl> wrote:
> First, your covariance matrix is not symmetric definite. I think you prob=
ably meant:
>
> CapSigma = {{sigma1^2, rho*sigma1*sigma2}, {rho*sigma1*sigma2,
>    sigma2^2}}
>
> Assuming that, you can get the conditional expectation from the definitio=
n:
>
> Integrate[y PDF[dist, {6.3, y}], {y, -Infinity, Infinity}]/
>  PDF[MarginalDistribution[dist, 1], 6.3]
>
> 5.6
>
> Alternatively, you can use Mathematica 8 built in NExpectation function:
>
> NExpectation[y \[Conditioned] 6.299 <= x <= 6.301, {x, y} \[Distribut=
ed] dist]
> 5.6
>
> Once you have the conditional expectation, you can compute the conditiona=
l variance, e.g.
>
>  Chop[Integrate[(y - 5.6)^2*PDF[dist, {6.3, y}], {y, -Infinity, Infinit=
y}]/PDF[MarginalDistribution[dist, 1], 6.3]]
>  0.0256
>
> Andrzej Kozlowski
>
> On 10 Jul 2011, at 11:01, Steve wrote:
>
>
>
> > Hello,
>
> > Can someone help me with this ?
>
> > I have 2 normal distributions; dist1 describes x and dist2 describes
> > y. Each are fully defined and are correlated to one another by the
> > correlation coefficient. How can I detemine the mean and standard
> > deviation of the expected normal distribution that is associated with
> > a given x value from dist1 ?
>
> > An example:
> > mean1 = 5.8
> > sigma1 =0 .2
>
> > mean2 = 5.3
> > sigma2 = 0.2
>
> > Correlation Coefficient, rho = 0.6
>
> > Given an x value of 6.3 (from dist1) what is the corresponding mean
> > and standard deviation of y ?
>
> > I can view the combined density function from the following:
>
> > Mu = {mean1, mean2}
> > CapSigma = {{sigma1^2, rho*sigma1*sigma2} , {rho,  rho*sigma1*sigma=
2}
> > dist = MultinormalDistribution[Mu,CapSigma]
> > pdf = PDF[dist,{x,y}]
> > plot1 = Plot3D[pdf, {x,4,7},{y,4,7}, PlotRange->All]
>
> > but can't see how to determine the mean and the standard deviation of
> > y for a given value of x, like 6.3
>
> > Any help would be appreciated.
>
> > Thanks,
>
> > --Steve- Hide quoted text -
>
> - Show quoted text -

Many thanks to all that responded to my question.

I apologize for the errors in my original post. I hadn't yet figured
out how to copy and paste Mathematica code into Usenet posts in a way
that is human readable, I therefore hand-typed my original post which
introduced typographical errors. Hopefully I can do better with this
and future posts.

The example problem of my original post came from a textbook on
mathematical statistics where their solution is mean = 5.6 and stddev
= 0.16

Andrzej Kozlowski's posting was particularly useful to me as it
confirmed the textbook solution and showed how I might generalize this
analysis to other situations. Below is my understanding of his
solution.

mean1 = 5.8;
sigma1 = .2;
mean2 = 5.3;
sigma2 = .2;
rho = .6;
mu = {mean1, mean2};
capsigma = {{sigma1^2, rho*sigma1*sigma2}, {rho*sigma1*sigma2,
sigma2^2}};
dist1 = MultinormalDistribution[mu, capsigma];
In[9]:= x1 = 6.3;
conditionalpdf = PDF[dist1, {x1, y}];
a = Integrate[y*conditionalpdf, {y, -Infinity, Infinity}]  ;
marginaldistribution = MarginalDistribution[dist1, 1];
b = PDF[marginaldistribution, x1] ;
mux2 = a/b;
c = Integrate[(y - mux2)^2*PDF[dist1, {x1, y}], {y, -Infinity,
    Infinity}];
sigmax2 = Sqrt[Chop[c/b]];
{mux2, sigmax2}
Out[18]= {5.6, 0.16}

What I really need to do is perform this analysis on test data for
which I have only a few data points, hence the Student T distribution
would be more appropriate than the Normal distribution. Secondly,
values for the "independent" and "dependent" variables have no
physical meaning below zero. So this implies that I need truncated
distributions. I'm hoping that the solution Andrzej  provided can be
generalized for these added complications.
Here are my 9 {F,t} data points where "F" is considered "independent"
and t considered "dependent".

{{1.01041, 0.3152}, {10.455, 0.3386}, {17.9032, 0.2534}, {24.9581,
   0.5412}, {26.4688, 0.3251}, {27.4651, 0.4428}, {30.1682,
   0.3402}, {36.6174, 0.2106}, {45.6129, 0.2154}}

Would someone be so kind as to plop this data into their notebook to
confirm a solution or two for me ? My results are below which are
based on truncating the Student T distribution, 8 degrees of freedom
and a calculated rho of -0.2327.

If F=1.0,    the probability of t exceeding 0.5 is 18.24%
If F=10.0,  the probability of t exceeding 0.5 is 13.17%
If F.0,   the probability of t exceeding 0.5 is 9.24%
If F=24.52,  the probability of t exceeding 0.5 is 8.11%
If F=30.0,   the probability of t exceeding 0.5 is 7.20%
If F=40.0,  the probability of t exceeding 0.5 is 6.55%

While working this, I came upon some integration problems and needed
to increase MaxRecursion to 20
Thanks in advance,
--Steve

[toc] | [prev] | [next] | [standalone]

#3803

On Jul 17, 3:03 am, Steve <s...@epix.net> wrote:
> [...]
> What I really need to do is perform this analysis on test data for
> which I have only a few data points, hence the Student T distribution
> would be more appropriate than the Normal distribution. Secondly,
> values for the "independent" and "dependent" variables have no
> physical meaning below zero. So this implies that I need truncated
> distributions. I'm hoping that the solution Andrzej  provided can be
> generalized for these added complications.
> Here are my 9 {F,t} data points where "F" is considered "independent"
> and t considered "dependent".
>
> {{1.01041, 0.3152}, {10.455, 0.3386}, {17.9032, 0.2534}, {24.9581,
>    0.5412}, {26.4688, 0.3251}, {27.4651, 0.4428}, {30.1682,
>    0.3402}, {36.6174, 0.2106}, {45.6129, 0.2154}}
>
> Would someone be so kind as to plop this data into their notebook to
> confirm a solution or two for me ? My results are below which are
> based on truncating the Student T distribution, 8 degrees of freedom
> and a calculated rho of -0.2327.
>
> [...]

I have several comments. First, the correlation of t with F is
so small that it is hard to justify treating it as nonzero. The
unbiased estimate of the conditional variance of t|F is bigger than
the unbiased estimate of the marginal variance of t. (The happens
whenever the F-statistic for testing the significance of the
correlation is < 1.)

In bivariate normal correlation, and in linear regression with
homoscedastic normal error, df = n-2, not n-1.

Regression models require only the conditional distribution of the
dependent variable given the independent variable. The independent
variable need not be random.

The fact that t can not be negative means that its conditional
distributions can not be normal. Is ordinary least squares fitting
justified? Yes, but only if conditional normality is abandoned. One
solution is to treat the conditional distributions as Gamma[a,b]
variables, where a is the shape constant and b is the scale constant.
Take a = m[F]^2/v and b = v/m[F]. Then the mean of each conditional
distribution will be m[F], the variance of each conditional
distribution will be v, and the Gauss-Markov theorem justifies
ordinary least-squares fitting.

Regardless of whether the conditional distributions are assumed
to be heteroscedastic truncated normal or homoscedastic gamma,
the sampling distribution of the estimates of the regression
coefficients and the conditional variance will not be the same
as in the usual homoscedastic normal case, and the usual Student-t
distributions can not be used to estimate quantiles of the conditional
distributions.

This may be a situation where one of John Tukey's antihubrisines
applies:
 "The data may not contain the answer. The combination of some
  data and an aching desire for an answer does not ensure that a
  reasonable answer can be extracted from a given body of data."

[toc] | [prev] | [next] | [standalone]

#3849

On Jul 17, 3:03 am, Steve <s...@epix.net> wrote:
> [...]
> What I really need to do is perform this analysis on test data for
> which I have only a few data points, hence the Student T distribution
> would be more appropriate than the Normal distribution. Secondly,
> values for the "independent" and "dependent" variables have no
> physical meaning below zero. So this implies that I need truncated
> distributions. I'm hoping that the solution Andrzej  provided can be
> generalized for these added complications.
> Here are my 9 {F,t} data points where "F" is considered "independent"
> and t considered "dependent".

> {{1.01041, 0.3152}, {10.455, 0.3386}, {17.9032, 0.2534}, {24.9581,
>    0.5412}, {26.4688, 0.3251}, {27.4651, 0.4428}, {30.1682,
>    0.3402}, {36.6174, 0.2106}, {45.6129, 0.2154}}

> Would someone be so kind as to plop this data into their notebook to
> confirm a solution or two for me ? My results are below which are
> based on truncating the Student T distribution, 8 degrees of freedom
> and a calculated rho of -0.2327.

> [...]

Another approach is to regress u = log[t] on f linearly. This
solves the problem of keeping the conditional distributions of
t non-negative, but makes the regression of t on f nonlinear.

For your data, switching to u increases the correlation, but
the conditional s.d. is still bigger than the marginal s.d.,
so there is still room for questioning the whole exercise.

FWTW, here are the numbers I got:

{mf, sf} = {Mean@f, StandardDeviation@f}

{24.5177, 13.3704}

{mu, su} = {Mean@u, StandardDeviation@u}

{-1.1482, .311145}

r = Correlation[f,u]

-.322776

b = r*su/sf  (* slope *)

-.00751142

a = mu - b*mf  (* intercept *)

-.964034

se = su*Sqrt[(1-r^2)(n-1)/(n-2)]  (* conditional s.d. *)

.314825

Table[{x, CDF[StudentTDistribution[n-2],
      (a+b*x-Log[.5])/se]}, {x,0,50,5}]
(* f, Pr[(t|f) > .5] *)

{{0, .209020},
{ 5, .179928},
{10, .154056},
{15, .131283},
{20, .111422},
{25, .0942435},
{30, .0794909},
{35, .0669001},
{40, .0562108},
{45, .0471758},
{50, .0395667}}

[toc] | [prev] | [standalone]

Back to top | Article view | comp.soft-sys.math.mathematica

csiph-web