Carga de funciones y librerías

# librerías
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   1.0.0 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.5.0 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(markovchain)
## Package:  markovchain
## Version:  0.9.0
## Date:     2022-07-01
## BugReport: https://github.com/spedygiorgio/markovchain/issues
# Matriz de probs. transición de n pasos
ptran.n=function(proceso,n){
  # proceso es La cmtd definida con markovchain
  # n son los pasos a dar
  nestados=length(states(proceso))
  #nestados=dim(proceso)
  i=1
  ptran=proceso[1:nestados]
  pn=ptran
  while(i<n){
    pn=pn%*%ptran
    i=i+1
  }
return(pn)  
}
#---------------------------------
# Cálculo de la matriz de ocupación en una CMTD
mocupa.proceso <- function(proceso, n)
{
  # proceso=proceso definido con markovchain
  # n son los pasos/transiciones a dar

  # Número de estados del proceso
  nestados <- dim(proceso)
  # Estados
  estados<- names(proceso)
  # Generamos la matriz de ocupaciones
  P <- proceso[1:nestados]
  # Inicialización con P^0=I
  mocupa=diag(nestados)
  for (i in 1:n)
   mocupa <- mocupa + ptran.n(proceso,i)
  
  return(mocupa)
}

Enunciado

Sabemos que en cada idioma las frecuencias de transición entre vocales y consonantes son diferentes. A partir de análisis recurrentes con textos de dos idiomas hemos identificado las probabilidades de transición entre vocales y consonantes en cada uno, estados=(vocal,consonante), y que vienen dadas a continuación.

\[p1=\left(\begin{matrix} 0.51 & 0.49 \\ 0.90 & 0.10 \end{matrix} \right), \qquad p2=\left(\begin{matrix} 0.25 & 0.75 \\ 0.30 & 0.70 \end{matrix} \right)\]

00. Descripción

Describe los procesos.

Definimos primero los procesos para la librería markovchain a partir de las probabilidades de transición.

estados=c("vocal","consonante")
# matriz de transición idioma1
p1=matrix(c(0.51,0.49,0.9,0.1),byrow = TRUE,ncol=2,dimnames=list(estados,estados))
# proceso 1: idioma1
idioma1=new("markovchain",states=colnames(p1),byrow=TRUE,transitionMatrix=p1,name="idioma1")
# matriz de transición idioma2
p2=matrix(c(0.25,0.75,0.3,0.7),byrow = TRUE,ncol=2,dimnames=list(estados,estados))
# proceso 2: idioma2
idioma2=new("markovchain",states=colnames(p2),byrow=TRUE,transitionMatrix=p2,name="idioma2")

# y lo mostramos en formato data.frame
as(idioma1,"data.frame")
##           t0         t1 prob
## 1      vocal      vocal 0.51
## 2      vocal consonante 0.49
## 3 consonante      vocal 0.90
## 4 consonante consonante 0.10
idioma1
## idioma1 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  vocal, consonante 
##  The transition matrix  (by rows)  is defined as follows: 
##            vocal consonante
## vocal       0.51       0.49
## consonante  0.90       0.10
as(idioma2,"data.frame")
##           t0         t1 prob
## 1      vocal      vocal 0.25
## 2      vocal consonante 0.75
## 3 consonante      vocal 0.30
## 4 consonante consonante 0.70
idioma2
## idioma2 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  vocal, consonante 
##  The transition matrix  (by rows)  is defined as follows: 
##            vocal consonante
## vocal       0.25       0.75
## consonante  0.30       0.70

Verificamos que todos sus estados son recurrentes y la cadena es irreducible.

# Describimos los procesos
summary(idioma1)
## idioma1  Markov chain that is composed by: 
## Closed classes: 
## vocal consonante 
## Recurrent classes: 
## {vocal,consonante}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE
summary(idioma2)
## idioma2  Markov chain that is composed by: 
## Closed classes: 
## vocal consonante 
## Recurrent classes: 
## {vocal,consonante}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

01. Transición de un paso

En cada idioma, ¿cuál es la probabilidad de encontrar dos vocales seguidas en un texto? ¿Cuál es la probabilidad de que si un texto empieza por vocal, el siguiente carácter sea consonante?

# A partir del proceso
cat("Respuesta idioma1:",idioma1[1:2][1,1])
## Respuesta idioma1: 0.51
# matriz de transición
p1[1,1]
## [1] 0.51
# probabilidades de transición
transitionProbability(idioma1,t0="vocal",t1="vocal")
## [1] 0.51
# Probabilidad condicional de un paso, dado un estado inicial
conditionalDistribution(idioma1,"vocal")
##      vocal consonante 
##       0.51       0.49
cat("Respuesta idioma2:",idioma2[1:2][1,1])
## Respuesta idioma2: 0.25

02. Transición de n pasos

En un texto de 10 caracteres, ¿qué probabilidad hay de que si empieza en vocal, acabe en consonante?

# Un texto de 10 caracteres tiene 9 transiciones
ptran.n(idioma1,9)
##                vocal consonante
## vocal      0.6474084  0.3525916
## consonante 0.6476172  0.3523828
# respuesta
cat("Idioma1:" , ptran.n(idioma1,9)[1,2])
## Idioma1: 0.3525916
ptran.n(idioma2,9)
##                vocal consonante
## vocal      0.2857143  0.7142857
## consonante 0.2857143  0.7142857
# respuesta
ptran.n(idioma2,9)[1,2]
## [1] 0.7142857
cat("Idioma2:" ,ptran.n(idioma2,9)[1,2])
## Idioma2: 0.7142857

03. Grafos

Representa sendos idiomas con un grafo.

plot(idioma1)

plot(idioma2)

04.Simulación y estimación MC

Por simulación aproxima (MC) el porcentaje de saltos vocal-consonante que se producen en un texto de 20 caracteres.

Algoritmo

Elegimos idioma. 1. Simular texto de 20 caracteres 2. ccontar el número de saltos de vocal a consonante y dividir por 19 (saltos) 3. Repetir pasos 1 y 2 nsim veces para conseguir nsim valores de la proporción de saltos 4. Calcular (MC) el promedio de la proporción saltos y el error asociado

nsim=5000
n_carac=20
p=c()
for(i in 1:nsim){
  # simulación secuencia de n_carac caracteres
  texto_1=markovchainSequence(n_carac,idioma1,include.t0 = FALSE)
  # contabilizamos saltos
  p[i]=createSequenceMatrix(texto_1)[2,1]/(n_carac-1)
}
estim=mean(p)
error=sd(p)/sqrt(nsim)
ic_low=estim-qnorm(0.975)*error
ic_up=estim+qnorm(0.975)*error
cat("Estimación idioma1: ",estim,"[" ,ic_low,",",ic_up,"]")
## Estimación idioma1:  0.3174316 [ 0.3156683 , 0.3191949 ]

05. Probabilidades marginales

¿Cuál es la probabilidad de que el tercer carácter en un texto sea una vocal, si sabemos que la proporción de palabras que empiezan por vocal en el idioma1 es de 0.5 y en el idioma2 de 0.2.

# distribución inicial
p1=c(0.5,0.5)
p2=c(0.2,0.8)

# Probabilidad marginal a 2 pasos
pmarginal1=p1%*%ptran.n(idioma1,2);pmarginal1
##        vocal consonante
## [1,] 0.62505    0.37495
pmarginal2=p2%*%ptran.n(idioma2,2);pmarginal2
##       vocal consonante
## [1,] 0.2855     0.7145
cat("\n Idioma1:", pmarginal1[1,1])
## 
##  Idioma1: 0.62505
cat("\n Idioma2:", pmarginal2[1,1])
## 
##  Idioma2: 0.2855

06.Tiempos de ocupación

En cada idioma, en una palabra de 5 caracteres que empieza por vocal, ¿cuántas vocales esperamos encontrar? ¿Y si la palabra empieza por consonante?

Nos pregunta por el número de visitas o tiempo de ocupación de cada uno de los estados (vocal/consonante) cuando llegamos a 5 caracteres (en 4 transiciones entre caracteres). Calculamos pues la matriz de tiempos de ocupación hasta el instante \(n=4\), para la que utilizamos la función mocupa.proceso().

n=4
mocupa1=mocupa.proceso(idioma1,n);mocupa1
##               vocal consonante
## vocal      3.493308   1.506692
## consonante 2.767393   2.232607
mocupa2=mocupa.proceso(idioma2,n);mocupa2
##               vocal consonante
## vocal      2.108844   2.891156
## consonante 1.156462   3.843537
cat("\n Idioma1 si empieza por vocal:",mocupa1[1,1])
## 
##  Idioma1 si empieza por vocal: 3.493308
cat("\n Idioma1 si empieza por consonante:",mocupa1[2,1])
## 
##  Idioma1 si empieza por consonante: 2.767393
cat("\n Idioma2 si empieza por vocal:",mocupa.proceso(idioma2,n)[1,1])
## 
##  Idioma2 si empieza por vocal: 2.108844
cat("\n Idioma2 si empieza por consonante:",mocupa.proceso(idioma2,n)[2,1])
## 
##  Idioma2 si empieza por consonante: 1.156462

07. Tiempos de primer paso

¿Cuántos caracteres habremos de leer por término medio en un texto (en cada idioma) hasta encontrar la primera vocal?

# Calculamos todos los tiempos de primer paso
meanFirstPassageTime(idioma1)
##               vocal consonante
## vocal      0.000000   2.040816
## consonante 1.111111   0.000000
meanFirstPassageTime(idioma2)
##               vocal consonante
## vocal      0.000000   1.333333
## consonante 3.333333   0.000000
# Puesto que se trata del primer paso, el primer carácter necesariamente ha de ser una consonante
cat("\n Respuesta Idioma1:",meanFirstPassageTime(idioma1)[2,1])
## 
##  Respuesta Idioma1: 1.111111
cat("\n Respuesta Idioma2:",meanFirstPassageTime(idioma2)[2,1])
## 
##  Respuesta Idioma2: 3.333333

08. Probabilidad de primer paso

En un texto que se inicia con una vocal, ¿con qué probabilidad encontraremos la primera consonante en el cuarto carácter? Compara los resultados para los dos idiomas.

# probabilidad de primer paso por el estado "consonante" partiendo de una "vocal"
# tras n=3 transiciones
n=3
fpm1=firstPassageMultiple(idioma1,state="vocal",set=c("consonante"),n=n);fpm1
##        set
## 1 0.490000
## 2 0.249900
## 3 0.127449
fpm2=firstPassageMultiple(idioma2,state="vocal",set=c("consonante"),n=n);fpm2
##        set
## 1 0.750000
## 2 0.187500
## 3 0.046875
cat("\n Resultado idioma1:",fpm1[3])
## 
##  Resultado idioma1: 0.127449
cat("\n Resultado idioma2:",fpm2[3])
## 
##  Resultado idioma2: 0.046875

09. Tiempo de recurrencia

A partir de una consonante, ¿cuántos caracteres hay que avanzar en el texto, por término medio, hasta encontrar otra consonante en cada idioma?

# tiempo medio de recurrencia, esto es, volver a encontrar otra consonante si partimos de una consonante
mrt1=meanRecurrenceTime(idioma1);mrt1
##      vocal consonante 
##   1.544444   2.836735
mrt2=meanRecurrenceTime(idioma2);mrt2
##      vocal consonante 
##        3.5        1.4
cat("\n Resultado idioma1:",mrt1[2])
## 
##  Resultado idioma1: 2.836735
cat("\n Resultado idioma2:",mrt2[2])
## 
##  Resultado idioma2: 1.4

10. Distribución estacionaria

¿Qué proporción de vocales y consonantes hay en cada idioma? En un texto de 1000 caracteres que empieza por vocal, ¿cuántas vocales esperamos encontrar? ¿Y si empieza por consonante?

Para contestar la primera pregunta recurrimos a la distribución estacionaria, que nos da la probabilidad estacionaria para cada uno de los estados posibles, o lo que es lo mismo, la proporción de vocales y consonantes. Para responder la segunda pregunta, puesto que estamos en un texto muy largo con 1000 caracteres, ya no utilizaremos a la matriz de ocupación hasta un instante \(n=1000\), sino la distribución de ocupación, que coincide con la distribución estacionaria.

Sabemos que una CMTD irreducible y aperiódica tiene una única distribución estacionaria, que coincide con la distribución de ocupación. Que es irreducible ya lo comprobamos anteriormente al definir el sistema; verifiquemos pues que es aperiódica, esto es, que su periodo es 1, con la función period().

period(idioma1)
## [1] 1
period(idioma2)
## [1] 1
#calculamos la distribución estacionaria
pi1=steadyStates(idioma1); pi1
##         vocal consonante
## [1,] 0.647482   0.352518
pi2=steadyStates(idioma2); pi2
##          vocal consonante
## [1,] 0.2857143  0.7142857

Tenemos que en el idioma1 el 65% de los caracteres en un texto son vocales, que predominan sobre las consonantes, mientras que en el idioma2 sólo un 29%, y la supremacía es de las consonantes, con un 71.4%.

Para calcular el número esperado de vocales y consonantes en un texto de 1000 caracteres, basta multiplicar estas probabilidades por 1000.

# número esperado de vocales y consonantes en un texto de 1000 caracteres
round(1000*pi1)
##      vocal consonante
## [1,]   647        353
round(1000*pi2)
##      vocal consonante
## [1,]   286        714

11. Ecuación de balance

Verifica que la distribución estacionaria verifica la ecuación de balance o del estado estacionario.

# es igual a su producto por la matriz de transición
pi1%*%idioma1[1:2]
##         vocal consonante
## [1,] 0.647482   0.352518
pi2%*%idioma2[1:2]
##          vocal consonante
## [1,] 0.2857143  0.7142857

12.Estado estacionario

Calcula secuencialmente la proporción de vocales que vamos a encontrar en un texto de 1,2,3, …,10, 20, 50, 100 caracteres, y relaciona los resultados con el estado estacionario. Utiliza la distribución inicial que dimos en el apartado 05.

# Nos piden la  distribución marginal a n pasos para n en la secuencia dada
# distribución inicial
p1=c(0.5,0.5)
p2=c(0.2,0.8)
secuencia=c(1:10,20,50,100)
for(i in secuencia){
# Probabilidad marginal a i pasos
  pmarginal1=p1%*%ptran.n(idioma1,i)
  cat("\n n=",i)
  cat("\n Idioma1:",pmarginal1[1])
  pmarginal2=p2%*%ptran.n(idioma2,i)
  cat("\n Idioma2:",pmarginal2[1])
}
## 
##  n= 1
##  Idioma1: 0.705
##  Idioma2: 0.29
##  n= 2
##  Idioma1: 0.62505
##  Idioma2: 0.2855
##  n= 3
##  Idioma1: 0.6562305
##  Idioma2: 0.285725
##  n= 4
##  Idioma1: 0.6440701
##  Idioma2: 0.2857137
##  n= 5
##  Idioma1: 0.6488127
##  Idioma2: 0.2857143
##  n= 6
##  Idioma1: 0.6469631
##  Idioma2: 0.2857143
##  n= 7
##  Idioma1: 0.6476844
##  Idioma2: 0.2857143
##  n= 8
##  Idioma1: 0.6474031
##  Idioma2: 0.2857143
##  n= 9
##  Idioma1: 0.6475128
##  Idioma2: 0.2857143
##  n= 10
##  Idioma1: 0.64747
##  Idioma2: 0.2857143
##  n= 20
##  Idioma1: 0.647482
##  Idioma2: 0.2857143
##  n= 50
##  Idioma1: 0.647482
##  Idioma2: 0.2857143
##  n= 100
##  Idioma1: 0.647482
##  Idioma2: 0.2857143

13. Ajuste de CMTD

Supón que accedes a un texto que transformas en una secuencia de vocales y consonantes (disponible para descarga en Github). Verifica, con dicha secuencia que proviene de una CMTD.

Para ello podemos utilizar la función verifyMarkovProperty() que resuelve un test de hipótesis basado en la Chi-cuadrado, sobre el cumplimiento de la propiedad de Markov con los datos proporcionados. Obtener un p-valor significativo significa que rechazaríamos esa propiedad y por lo tanto rechazaríamos que se trata de una CMTD.

# leemos los datos de Github
texto=read.csv("https://raw.githubusercontent.com/UMH1477/data1477/main/textos.csv")
# visualizamos el formato de los datos, con 2 columnas de texto
head(texto)
##       texto1     texto2
## 1      vocal consonante
## 2      vocal consonante
## 3 consonante consonante
## 4      vocal      vocal
## 5 consonante consonante
## 6 consonante consonante
# y ejecutamos el test sobre una de las columnas de texto
verifyMarkovProperty(texto$texto1)
## Testing markovianity property on given data sequence
## Chi - square statistic is: 2.056533 
## Degrees of freedom are: 5 
## And corresponding p-value is: 0.8412687

El p-valor es de 0.84, que no permite rechazar la propiedad de Markov.

Por último, para testar la estacionariedad, esto es, si \(p_{ij}(n)=p_{ij}\) para cualquier \(n\) suficientemente grande, tenemos la función assessStationarity(). Obtener un p-valor significativo significa que rechazaríamos esa propiedad y por lo tanto rechazaríamos que se trata de una CMTD.

14. Contraste de similaridad

Con los datos utilizados en el apartado anterior, verifica si alguno de los textos que se proporcionan en las columnas de la base de datos, se podría asimilar como perteneciente a alguno de los idiomas presentados, idioma1 o idioma2.

Para contrastar la similitud entre una secuencia y un proceso teórico que tenemos definido podemos utilizar la función verifyEmpiricalToTheoretical(). Veamos cómo funciona contrastando si cualquiera de las secuencias en la bd anterior es compatible con el idioma1 y el idioma2.

Para que funcione adecuadamente, hemos de convertir a numéricas las etiquetas de los estados.

states(idioma1)
## [1] "vocal"      "consonante"
# estados en orden "vocal"-"consonante"
# asignamos 0 a vocal, y 1 a consonante
p_idioma1=idioma1[1:2]
dimnames(p_idioma1)=list(0:1,0:1)

texto1=(texto$texto1=="consonante")*1
texto2=(texto$texto2=="consonante")*1

# compatibilidad con idioma1
verifyEmpiricalToTheoretical(texto1,as(p_idioma1,"markovchain"))
## Testing whether the
##       0     1
## 0 32800 31854
## 1 31854  3491
## transition matrix is compatible with
##      0    1
## 0 0.51 0.49
## 1 0.90 0.10
## [1] "theoretical transition matrix"
## ChiSq statistic is 2.460873 d.o.f are 2 corresponding p-value is 0.292165
## $statistic
##        0 
## 2.460873 
## 
## $dof
## [1] 2
## 
## $pvalue
##        0 
## 0.292165
# texto1: no rechaza compatibilidad con idioma1
verifyEmpiricalToTheoretical(texto2,as(p_idioma1,"markovchain"))
## Testing whether the
##       0     1
## 0  7198 21599
## 1 21599 49603
## transition matrix is compatible with
##      0    1
## 0 0.51 0.49
## 1 0.90 0.10
## [1] "theoretical transition matrix"
## ChiSq statistic is 153716.3 d.o.f are 2 corresponding p-value is 0
## $statistic
##        0 
## 153716.3 
## 
## $dof
## [1] 2
## 
## $pvalue
## 0 
## 0
# texto2: rechaza compatibilidad con idioma1 a favor de diferencias

15. Análisis de costes

En un concierto de música, el organizador decide vender las entradas a dos precios: 50€ y 70€. Para ello, etiqueta los tiquets que va vendiendo, con una vocal o una consonante, en función de una secuencia de vocales-consonantes que ha simulado utilizando la CMTD del idioma1. Los tiquets con etiqueta “ vocal” los cobra a 50€ y los etiquetados con “ consonante” a 70€.

1. ¿Cuánto dinero habrá ingresado por la venta de tiquets cuando haya vendido 20 entradas?

benef=c(50,70)
# Beneficios a horizonte finito
mocupa.proceso(idioma1,19)%*%benef
##                [,1]
## vocal      1135.935
## consonante 1150.323
sol=round(expectedRewards(idioma1,19,benef),2)

Si la primera entrada que vende es de 50€, ingresa 1135.93€; si es de 70€, ingresa 1150.32€.

1. ¿Cuánto dinero habrá ingresado por término medio por entrada vendida cuando cierre caja (se espera que venda más de 5000 entradas)? ¿Y cuando haya vendido 5000 entradas?

# Beneficio medio a largo plazo
ganancia_media=round(steadyStates(idioma1)%*%benef,2)
ganancia_5000=round(steadyStates(idioma1)%*%benef*5000,2)

Habrá ganado de media unos 57.05€ por entrada. Con 5000 entradas habrá ingresado un total de 2.852518^{5}€.