| #- EST0116 INFERENCIA II | |
| #- Slides 0 | |
| # slide 23 ---------------------------------------------------------------- | |
| estu_sim1 <- function(eps, theta, K = 100) { | |
| ns <- seq(5, 10000, 1000) | |
| prob <- numeric() | |
| for (n in ns) { | |
| cont <- numeric() |
| library(tidyverse) | |
| library(ggplot2) | |
| diamonds %>% | |
| count(cut) %>% | |
| arrange(desc(cut)) %>% | |
| mutate(prop = n/sum(n) ) %>% | |
| mutate(lab.ypos = cumsum(prop) - 0.5*prop) %>% | |
| ggplot(aes(x = "", y = prop, fill = cut)) + | |
| geom_bar(width = 1, stat = "identity", color = "white") + |
| rNx <- function(k = 1E3, x = 1) { | |
| out <- numeric() | |
| for (i in 1:k) { | |
| nx <- 1 | |
| u <- runif(1) | |
| while(u <= x) { | |
| u <- u + runif(1) | |
| nx <- nx + 1 | |
| } | |
| out <- append(out, nx) |
| # Adaptacao da funcao 'mwar.ani' do pacote 'animation' do R para | |
| # ilustrar o conceito de media movel para k observacoes consecutivas | |
| mm.ani = function(x, k = 7, mat = matrix(1:2, 2), | |
| widths = rep(1, ncol(mat)), | |
| heights = rep(1, nrow(mat)), lty.rect = 2, ...) { | |
| nmax = ani.options('nmax') | |
| n = length(x) | |
| if (k > n) | |
| stop('A largura da janela k deve ser menor que o tamanho de x!') |
These are the R codes to implement the illustrative example in Costa, Paulino & Singer (2020, submitted).
The MIT License (MIT)
Copyright (c) 2020 Eliardo G. Costa, Carlos Daniel Paulino & Julio M. Singer
Permission is hereby granted, free of charge, to any person obtaining a copy of
| n <- 50 | |
| w <- 1 | |
| set.seed(161017) | |
| lam <- sample(c(5, 25, 50), n, replace = TRUE, prob = c(0.35, 0.3, 0.35)) # 3 clusters! | |
| x <- rpois(n, w*lam) | |
| hist(x, prob = TRUE, breaks = 15) | |
| source("https://gist.githubusercontent.com/eliardocosta/7c350eece9f09a7caa0f48852621063c/raw/dd6eae6d9bb3101bdeeff0d9d087db91f9d7e74d/ppostPoiDP.R") | |
| foo <- ppostPoiDP(x = x, w = w, phi = sqrt(mean(x)/var(x)), lam0 = mean(x), alpha = 1, nburn = 500, nsam = 1E3) | |
| hist(foo$nclu) # histograma num. de clusters por iteracao | |
| plot(foo$lam, foo$fdist) # distribuicao acumulada a posteriori |
| r <- 0.5 # raio da circunferencia | |
| n <- 1E5 | |
| x <- runif(n) | |
| y <- runif(n) | |
| pdentro <- numeric(n) | |
| for (i in seq_len(n)) { # contando os pontos dentro da circunferencia | |
| pdentro[i] <- ifelse((x[i] - 0.5)^2 + (y[i] - 0.5)^2 <= r^2, 1, 0) | |
| } |
| {fileheader} | |
| main <- function() { | |
| return(cat()) | |
| } | |
| if(interactive()) main() |
1- Go to Tools -> Global Options -> Code -> Snippets -> Edit Snippets.
2- Then paste the following code:
snippet 0bsd
`r paste('# script.R')`
`r paste('#')`
`r paste('# Copyright (c)', format(Sys.time(), "%Y"), '<Your name>')`
`r paste('#')`
`r paste('# This program is free software; you can redistribute it and/or modify')`
`r paste('# it under the terms of the BSD Zero Clause License as published by')`
1- Go to Tools -> Global Options -> Code -> Snippets -> Edit Snippets.
2- Then paste the following code:
snippet mit
`r paste('# script.R')`
`r paste('#')`
`r paste('# The MIT License (MIT)')`
`r paste('#')`
`r paste('# Copyright (c)', format(Sys.time(), "%Y"), '<Your name>')`
`r paste('#')`