1 人文2019問1

1.1 条件付き期待値 \(E(Z|Z>0)\).

　　　\(\displaystyle E(Z|Z\geq0) =\int ^{\infty }_{0}z\frac{\varphi\left(z\right)}{P(Z\geq 0)}dz =2\int ^{\infty }_{0}z\varphi\left(z\right)dz\).
　　\(\therefore E(Z|Z\geq0)=2\left[-\varphi \left( z\right) \right]_{0}^{\infty} =2\varphi(0)=2\frac{1}{\sqrt{2\pi}}=\sqrt{\frac{2}{\pi}} \fallingdotseq 0.79\).
乱数から条件を満たす要素を選び、平均をとる。

X <- rnorm(100)
head(X[X>0])

## [1] 0.5372363 1.6625520 1.1050737 0.3348966 0.4730636 0.9710500

mean(X[X>0])#条件付き期待値

## [1] 0.9348875

sqrt(2/pi)#理論値

## [1] 0.7978846

1.2 \(E(Z^2|Z>0)\).

\(\begin{split} E(Z^2|Z\geq0) &=\int ^{\infty }_{0}z^2\frac{\varphi\left(z\right)}{P(Z\geq 0)}dz =2\int ^{\infty }_{0}z^2\varphi\left(z\right)dz\\ &=2\int ^{\infty }_{0}z\cdot z\varphi\left(z\right)dz =2[z(-\varphi(z)]_{0}^{\infty}+2\int_{0}^{\infty}\varphi(z)dz=1 \end{split}\).
　　　

X <- rnorm(10000)
Y <- X[X>0]^2
mean(Y)#条件付き期待値

## [1] 1.000419

理論値である 1 に近い値である。

1.3 条件付き分散 \(V(Z|Z>0)\).

\(V(Z|Z\geq 0)=1-\frac{2}{\pi}\).

X <- rnorm(1000)
var(X[X>0])

## [1] 0.335238

1-2/pi

## [1] 0.3633802

2 人文2016問2

条件付き期待値も線型性を持つ。標準正規分布に直して計算する。
\(X=\sigma Z+\mu\) 　（\(X\geq a\) のとき \(Z\geq c\) ）とすると、
　　　\(\begin{split} E(X|X\geq a)&=E(\sigma Z+\mu|Z\geq c)\\&=\sigma E(Z|Z\geq c)+E(\mu|Z\geq c)\\ &=\sigma \dfrac{\varphi(c)}{P(Z\geq c)}+\mu \end{split}\).
が一般に成立するので、

　　　\(E(X|X\geq 500) =80\dfrac{\varphi(0.4)}{P(Z\geq 0.4)}+468 =80\dfrac{\varphi(0.4)}{0.3446}+468\).
\(\varphi \left( 0.4\right) =\dfrac{1}{\sqrt{2\pi }}\exp(-0.08)=\dfrac{1}{\sqrt{2\pi}\cdot 1.0833}\)　であるから、

　　　\(E(X|X\geq 500)=\dfrac{80}{\sqrt{2\pi}(1.0833)(0.3446)}+468=\dfrac{80}{0.9355}+468=553.5\)

X <- rnorm(10000,468,80)
head(X[X>500])

## [1] 588.8436 546.0415 591.3025 535.4339 614.6297 506.8617

mean(X[X>500])#条件付き期待値

## [1] 553.724

乱数による結果は理論値と近いものとなっている。

3 人文2017問3

library(psych)
library(MASS)
mu <- c(50,50)
Sigma <- matrix(c(15^2, 135,135,15^2), 2, 2)
data <- mvrnorm(100, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2])
attach(d)

## The following objects are masked _by_ .GlobalEnv:
## 
##     X, Y

head(d)

##          X        Y
## 1 27.21986 59.93850
## 2 32.04736 46.93937
## 3 38.06400 45.68251
## 4 58.98686 69.14100
## 5 29.72472 40.44529
## 6 50.41184 66.16211

(4).
\(E(Y|X\geq 50)=\alpha+\beta\cdot E(X|X\geq 50)\)

dc <- d[d$X>=50, ]
mean(dc$X)

## [1] 60.52917

mean(dc$Y)

## [1] 58.89275

https://multivariate-statistics.com/2022/06/09/r-programming-correlation-coefficient-plot/

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(scatterplot3d)
library(mvtnorm)
library(tidyr)
library(gapminder)
options(rgl.printRglwidget = TRUE) 
library(rgl)

scatterplot3d(d[,1], d[,2], dmvnorm(d, mean=c(50,50), sigma=Sigma),highlight=TRUE)

x <- d$X
y <- d$Y
z <- dmvnorm(d, mean=c(50,50),sigma=Sigma)

open3d()

## glX 
##   1

plot3d(x, y, z, type = "s", col = "blue", size = 1)

http://leeswijzer.org/R/R-binormal.html

library(mvtnorm)
n <- 20
mu <- c(0, 0)
rhos <- c(0.01, 0.5, 0.8, 0.99)
 
mkMatrix <- function(rho) {
  return(matrix(c(1, rho, rho, 1), ncol = 2))
}
 
#確率密度
par(mfrow = c(2, 2))
par(mar = c(4, 3, 1, 1))
par(oma = c(0, 0, 0, 0))
par(mgp = c(2, 1, 0))
for (rho in rhos) {
  Sigma <- mkMatrix(rho)
  x <- seq(-3, 3, 0.1)
  y <- x
  f <- function(u, v) {
    c(dmvnorm(matrix(c(u, v), ncol = 2), mu, Sigma))
  }
  density <- outer(x, y, f)
  persp(x, y, density, theta = 0, phi = 60, expand = 0.5,
        xlim = c(-3, 3), ylim = c(-3, 3), sub =  paste0("ρ=", rho))
}