Download Statistical Analysis of Binomial, Gamma, and Poisson Distributions and more Assignments Statistics in PDF only on Docsity!
ST 762, HOMEWORK 3 EXTRA PROBLEMS SOLUTIONS, FALL 2007
- By a conditioning argument, the cdf of X is F (x) = Φ 1 (x)(1 − α) + Φb(x)α, where
Φb(x) =
∫ (^) x
−∞
2 πb^2
e−u
(^2) /(2b (^2) ) du,
where the integrand is the normal density φb(x) with mean zero and variance b^2. Thus, the density is f (x) = φ 1 (x)(1 − α) + φb(x)α. It is straightforward to see that E(X) = 0, E(X^2 ) = (1 − α) + αb^2 , E(X^3 ) = 0, and
E(X^4 ) = (1 − α)
∫ (^) ∞
−∞
x^4 φ 1 (x) dx + α
∫ (^) ∞
−∞
x^4 φb(x) dx,
which, by a change of variables in the second integral and using Problem 2(a), reduces to E(X^4 ) = 3(1 − α) + 3b^4 α = 3{(1 − α) + αb^4 }. Thus, for ǫ = X{(1 − α) + αb^2 }−^1 /^2 , E(ǫ) = 0, E(ǫ^2 ) = 1, E(ǫ^3 ) = 0, and
var(ǫ^2 ) = E(ǫ^4 ) − 1 =
3 {(1 − α) + αb^4 } {(1 − α) + αb^2 } − 1 = 2 + κ,
so that κ = 3
{ 1 − α + αb^4 (1 − α + αb^2 )^2
} .
(b) Substitution with α = 0.01 and α = 0.05 gives the following values of κ
b α = 0. 01 α = 0. 05 2 0.252 0. 3 1.630 4. 4 5.053 10.
Thus, for b = 3, α = 0.01, for instance, var(ǫ^2 ) = 3.630 vs. 2 for the normal. The excess kur- tosis increases fairly quickly. Even for a very small proportion of outliers like 0.01, the excess kurtosis can depart considerably from that of the normal. This suggests that deviation from normality caused by outliers may be important when estimating β by quadratic estimating equations, as we will see later.
- (a) One way to do this is to use the moment generating function of the standard normal, given by m(t) = et (^2) / 2 . Taking derivatives gives m′(t) = tet (^2) / 2 , m′′(t) = et (^2) / 2 + t^2 et (^2) / 2 , m′′′(t) = 3tet^2 /^2 + t^3 et^2 /^2 , and miv^ = 3et^2 /^2 + 6t^2 et^2 /^2 + t^4 et^2 /^2. Thus, E(ǫ^3 ) = m′′′(0) = 0 and E(ǫ^4 ) = miv(0 = 3, so that var(ǫ^2 ) = E(ǫ^4 ) − {E(ǫ^2 )}^2 = 3 − 1 = 2, and κ = 0. (b) The gamma is a scaled exponential family random variable, so we can use the results in Homework 2, Problem 1, as E(ǫ^3 ) = E{(Y − μ)^3 /(σ^2 μ^2 )^3 /^2 } = m 3 /m^3 /^2 = ζ and var(ǫ^2 ) = E(ǫ^4 ) − 1 = E{(Y − μ)^4 /(σ^2 μ^2 )^2 } − 1 = m 4 /m^22 − 1 , so that κ = m 4 /m^22 − 3. We wish to find these values for the gamma. Now from Homework 2 Extra Problems, Problem 1, recall that
ζ = σbξξξ(ξ)/{bξξ(ξ)}^3 /^2 ,
κ = σ^2 bξξξξ(ξ)/{bξξ (ξ)}^2.
For the gamma, b(ξ) = − log(−ξ), bξ (ξ) = − 1 /ξ, bξξ(ξ) = 1/xi^2 , bξξξ(ξ) = − 2 /ξ^3 , and bξξξξ(ξ) = 6/ξ^4. As ξ = − 1 /μ, we obtain ζ = 2σ and κ = 6σ^2. Note that the skewness and excess kurtosis are of orders σ and σ^2 , a fact that will be of interest to us later. (c) We are interested in finding the coefficient of skewness ζ and the excess kurtosis κ for a binomial random variable with mean kp and variance kp(1−p). We can use the same approach as in (b) here, as this is a scaled exponential family; in this case, it is easy to show that b(ξ) = k log(1+eξ), so that bξ(ξ) = keξ^ /(1+eξ) = kp, bξξ(ξ) = keξ^ /(1+eξ){ 1 −eξ/(1+eξ)} = kp(1−p),
bξξξ(ξ) = keξ 1 + eξ^
3 ke^2 ξ (1 + eξ^ )^2
2 ke^3 ξ (1 + eξ^ )^3 = kp(1 − p)(1 − 2 p).
bξξξξ(ξ) = keξ 1 + eξ^
7 ke^2 ξ (1 + eξ^ )^2
12 ke^3 ξ (1 + eξ^ )^3
6 ke^4 ξ (1 + eξ^ )^4
= kp(1 − p)(1 − 6 p + 6p^2 ).
Thus, substituting into the expressions for ζ and κ, we obtain
ζ = 1 − 2 p {kp(1 − p)}^1 /^2
, κ = 1 − 6 p + 6p^2 kp(1 − p)
Note here that the skewness and excess kurtosis in fact depend on the mean kp. This shows that, for regression models and this distribution, the third and fourth moments of ǫj |xj must depend on the mean; hence, the assumption that the ǫj are independent of the xj would be inappropriate.
- (a) If we define as in the notes
sj =
( Yj (Yj − fj )^2
) , mj =
( fj σ^2 g j^2
) ,
using obvious shorthand notation, then we have
Dj =
( fβj 2 σ^2 g^2 j νβj 0 2 σg^2 j
)
as in the notes. However, now we have
V (^) j =
( σ^2 g^2 j ζj σ^3 g j^3 ζj σ^3 g j^3 (2 + κj )σ^4 g^4 j
) .
Now V − j 1 =
|V (^) j |
( (2 + κj )σ^4 g^4 j −ζj σ^3 g^3 j −ζj σ^3 g^3 j σ^2 g j^2
) ,
|V (^) j | = σ^6 g^6 j (2 + κj − ζ^2 j ).
Thus, the equation may be written as ∑^ n
j=
( fβj 2 σ^2 g^2 j νβj 0 2 σg^2 j
) 1 σ^6 g^6 j (2 + κj − ζ j^2 )
( (2 + κj )σ^4 g^4 j −ζj σ^3 g j^3 −ζj σ^3 g^3 j σ^2 g^2 j
) ( Yj − fj (Yj − fj )^2 − σ^2 g^2 j
) = 0.
Multiplying out and simplifying, we eventually arrive at the joint estimating equations ∑^ n
j=
(2 + κj )fβj /(σgj ) − 2 νβj ζj σgj (2 + κj − ζ j^2 ) (Yj − fj ) +
2 νβj − ζj fβj /(σgj ) σ^2 g^2 j (2 + κj − ζ^2 j ) {(Yj − fj )^2 − σ^2 g j^2 } = 0,
(d) Under these conditions, we have σ = 0.5, gj = fj , ζj = 1, and κj = 1.5. Then νβj = fβj /fj , so that aj and cj simplify to
aj =
( 4 fβj /f (^) j^2 − 16 / 5
) , cj =
( 0 32 / 5 j
) .
The first equation, which corresponds to β, thus reduces to
∑^ n
j=
f (^) j− 1 (Yj − fj )fβj = 0 ,
where we may ignore the constant “4,” where this equation is linear in Yj. (e) We see that, in all these cases, the quadratic estimation equation for β reduces to the usual linear maximum likelihood equations, which have the GLS form. The binomial (b), gamma (c), and Poisson (d) distributions are members of the scaled exponential family class, so we know that the maximum likelihood estimating equations for β must have this property. For distributions in this class, the variance is completely determined by the mean. Thus, an obvious explanation for the way that this turned out is that trying to gain information on β from the variance via a quadratic equation is pointless when the first four moments of the data, which are used in forming the quadratic equation, are the same as those of a distribution in the scaled exponential family class (the binomial here). We can’t improve on the maximum likelihood estimator even if we try to do this. We will discuss this more later in the course.
- Expanding (1) about α∗^ “close to” α gives
∑^ n
j=
DTj (α∗)V − j 1 (α∗){sj (α∗) − mj (α∗)}
∑^ n
j=
∂/∂α {DTj (α∗)V − j 1 (α∗)}{sj (α∗) − mj (α∗)}
∑^ n
j=
DTj (α∗)V − j 1 (α∗)∂/∂α {sj (α∗) − mj (α∗)}
(^) (α − α∗).
Now ∂/∂α sj (α) = ∂/∂αsj [{Yj − f (xj , β)}^2 ], where this notation is meant to indicate that sj (α) is a function of α only through this argument. By the chain rule, letting “′” denote differentiation with respect to the argument, we have
s′ j [{Yj − f (xj , β)}^2 ] ∂/∂α, [{Yj − f (xj , β)}^2 ]
= − 2 {Yj − f (xj , β∗)}∂/∂z sj [{Yj − f (xj , β∗)}^2 ]∂/∂α {f (xj , β∗)}. Substituting this into the above and using the definition of Dj (α), we obtain
∑^ n
j=
DTj (α∗)V − j 1 (α∗){sj (α∗) − mj (α∗)}
∑^ n
j=
∂/∂α {DTj (α∗)V − j 1 (α∗)}{sj (α∗) − mj (α∗)}(α − α∗)
∑^ n
j=
DTj (α∗)V − j 1 (α∗)Dj (α∗)(α − α∗)
∑^ n
j=
{Yj − f (xj , β∗)}DTj (α∗)V − j 1 (α∗)s′ j (α∗)∂/∂α {f (xj , β∗)}(α − α∗).
The term on the second line is negligible by virtue of the fact that E(sj (α∗)|xj ) ≈ mj (α∗), so that this term contains the product of two “small” terms. The term in the fourth line depends on the product (α − α∗){Yj − f (xj , β∗)}, which is also “small.” Disregarding these terms yields the approximation
∑^ n
j=
DTj (α∗)V − j 1 (α∗){sj (α∗) − mj (α∗)}
∑^ n
j=
DTj (α∗)V − j 1 (α∗)Dj (α∗)(α − α∗),
which may be rearranged to give (6.14). The rest of the problem follows by the same manip- ulations as in the notes.
- (a) First, we need to determine γ. Note that γ must satisfy
E(ǫ^2 ) =
2 γ
∫ (^) ∞
−∞
ǫ^2 e−|ǫ|/γ^ dǫ
γ
∫ (^) ∞
0
ǫ^2 e−ǫ/γ^ dǫ
= 2 γ^2
either by brute-force integration by parts or by noting that this is the second moment of an exponential distribution. To ensure E(ǫ^2 ) = 1, then, we obtain that γ = 1/
We may now write down the form of the density of Yj |xj ; the joint density will be the product. We have for this choice of γ f (ǫ) =
√^1
e−
√ 2 ǫ ,
and ǫ = (Y − f )/(σg) in shorthand, so that the density of Y (given x, suppressed here) has the form fY (y) = f {(Y − f )/(σg)}
σg
Thus, the joint density of interest is
∏^ n
j=
√^1
2 σg(β, θ, xj )
exp
{ −
|Yj − f (xj , β)| σg(β, θ, xj )
} .
(b) First, we need E(|ǫ|) for the double exponential. From (a), this is, with γ = 1/
E(|ǫ|) =
2 γ
∫ (^) ∞
−∞
|ǫ|e−|ǫ|/γ^ dǫ =
∫ (^) ∞
0
ue−
√ 2 udu = 1/√ 2.
To obtain the “trick,” note from above that η must satisfy
eλη^ = n−^1
∑^ n j=
|Yj − f (xj , β)|λ gλ(β, θ, xj )
Substituting this into log L in (1) gives the profile loglikelihood
log Lmax = −(n/λ) log
n
− 1 ∑^ n
j=
|Yj − f (β, xj )|λ gλ(β, θ, xj )
−^ (n/λ)
∑^ n
j=
log gλ/n(β, θ, xj ) − n/λ.
Defining ˙g(β, θ, xj ) as in the notes, we may write
log Lmax = −(n/λ) log
n−^1
∑^ n j=
|Yj − f (β, xj )|λ^ g˙λ(β, θ, xj ) gλ(β, θ, xj )
− n/λ + (n/λ) log n,
so that maximizing log L in θ is equivalent to minimizing
∑^ n
j=
{ |Yj − f (β, xj )|λ/^2 g˙λ/^2 (β, θ, xj ) gλ/^2 (β, θ, xj )
} 2 .
Thus, we may use nonlinear regression procedures to maximize log L, regressing “dummy” data equal to zero for all j on
Fj (θ) =
|Yj − f (β, xj )|λ/^2 g˙λ/^2 (β, θ, xj ) gλ/^2 (β, θ, xj )
For the power variance model g(β, θ, xj ) = f θ(xj , β), this becomes
Fj (θ) = |Yj − f (β, xj )|λ/^2
{ f˙ (β) f (xj , β)
}θλ/ 2 .
- It is straightforward to see that differentiation of (7.11) with respect to σ yields
σ^2 = (n − p)−^1
∑^ n
j=
{Yj − f (xj , βˆ)}^2 g^2 (βˆ, θ, xj )
Substituting this in (7.11) gives the profile objective function
(n − p) 2
(n − p) 2
log
(n − p)−^1
∑^ n
j=
{Yj − f (xj , βˆ)}^2 g^2 (βˆ, θ, xj )
(^) − (n^ −^ p) 2
log{ g˙(θ)}(2n)/(n−p)}
(n − p) 2 log |N (βˆ, θ)|^1 /(n−p),
where ˙g(θ) =
∏n j=1 g 1 / (^2) (βˆ, θ, xj ). Note then that θ maximizing this profile objective function must maximize
(n − p) 2
log
∑^ n j=
{Yj − f (xj , βˆ)}^2 g^2 (βˆ, θ, xj )
g˙(θ)}(2n)/(n−p)|N (βˆ, θ)|^2 /{2(n−p)}
(^).
Ignoring constants, then, the REML estimator for θ must minimize
∑^ n
j=
[ {Yj − f (xj , βˆ)} g(βˆ, θ, xj )
g˙(θ)n/(n−p)|N (βˆ, θ)|^1 /2(n−p)
] 2 .
Thus
Fj (θ) = {Yj − f (xj , βˆ)} g(βˆ, θ, xj )
g˙(θ)n/(n−p)|N (βˆ, θ)|^1 /2(n−p).
