ito4303

library(ggplot2)
library(purrr)
library(dplyr)
library(lawstat)

# Normal distribution
test_norm <- function(R = 2000, N = c(30, 30),
                      mu = c(0, 0), sigma = c(1, 1)) {
  y <- purrr::map(seq_len(R), function(i) {
    x1 <- rnorm(N[1], mean = mu[1], sd = sigma[1])

ito4303

---
title: "Stein estimator"
author: "ITÔ, Hiroki"
date: "2018/3/22"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

ito4303

#!/usr/bin/env ruby -Ku
#%%%{CotEditorXInput=AllText}%%%
#%%%{CotEditorXOutput=ReplaceAllText}%%%

require 'cgi'

while $stdin.gets
  print CGI.escapeHTML($_)
end

ito4303

library(Rcpp)
library(inline)

## C++ source
mandelbrot.src <- "
Rcpp::NumericVector x(xx);
Rcpp::NumericVector y(yy);
unsigned short t = Rcpp::as<unsigned short>(threshold);
int nx = x.size(), ny = y.size();
Rcpp::NumericVector u(nx);

ito4303

data {
  int<lower = 0> N_ind;
  int<lower = 0> N_t;
  int<lower = 0, upper = 2> Y[N_ind, N_t]; // 0: missing value
}

parameters {
  real<lower = 0, upper = 1> phi;
  real<lower = 0, upper = 1> psi;
  real<lower = 0, upper = 1> p;

ito4303

library(purrr)

sim_inc_sample <- function(N_start = 20, N_end = 100,
                           mu = 0, sigma = 1, sig = 0.05) {
  not_sig <- TRUE
  N <- N_start
  x1 <- rnorm(N_end, mu, sigma)
  x2 <- rnorm(N_end, mu, sigma)
  while(not_sig & N <= N_end) {
    p_value <- t.test(x1[1:N], x2[1:N])$p.value

ito4303

## ANOVA and Kruskal-Wallis test for data with unequal variances

library(ggplot2)
library(purrr)
library(dplyr)

# Normal distribution
test_norm <- function(R = 2000, N = c(30, 30, 30),
                      mu = c(0, 0, 0), sigma = c(1, 1, 1)) {
  y <- purrr::map(seq_len(R), function(i) {

ito4303

# see https://qiita.com/th1nkd0g/items/3182e210b6a8b0ba0e79
env PYTHON_CONFIGURE_OPTS="--enable-framework CC=clang" pyenv install 3.6.5
pyenv global 3.6.5
sudo ln -s /usr/local/var/pyenv/versions/3.6.5/Python.framework/Versions/3.6 /Library/Frameworks/Python.framework/Versions/

ito4303

library(nimble)
library(coda)

## Data
k <- 10
x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(1, 2, 2, 6, 4, 5, 8, 9, 9, 9, 10)
n <- length(x)

## Code

ito4303

library(nimble)
library(dplyr)

# Data
set.seed(123)
n <- 40
omega <- 0.4
lambda <- 3
y <- rbinom(n, 1, omega) * rpois(n, lambda)