[貝氏] 煩死人的 Variational Bayes 初探

先說，搞懂這個東西的過程令人崩潰的程度大概是這樣

羽生結弦，日本天才花式滑冰選手
圖片來源http://abcnews.go.com/Sports/photos/tough-straight-face-midair-22463006/image-22525718

貝氏統計推論中，有些Posterior實在太複雜太難算，用sampling (如MCMC)抽一大堆貌似合理的樣本來進行參數推論是經常使用的替代方式。然而sampling也有既存的問題，比如說常常讓運算速度不夠快，或是完全來自隨機以致收斂性等表現難以掌握等等。總之，學界開始尋找其他逼近(approximation)的方法在某些場合下取代MCMC。

Variational Bayes算是最紅的，許多貝氏建模的文章都會提供相對應的演算細節。

大致上的概念是這樣的，我們感興趣的東西是$p\left ( \theta \right |x)$，並且

$p\left ( \theta \right |x) \propto p\left ( \theta \right )p\left ( x |\theta \right)$

找一個替代品$q\left(\theta\right)$，並且讓這個替代品和posterior很像很像，換句話說就是努力把

$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$

變得很小很小(沒錯，Kullback-leibler divergence這個小壞壞又出現了)

然後，因為我們有好多個parameter要推論，mean field assumption會是一個把情況變得簡單的假設，也就是所有的參數都是獨立的。也就是

$q\left(\theta \right )=\prod_iq\left(\theta_i \right )$

把這個該變小的Kullback-leibler divergence寫開

$KL \left[q \left( \theta \right),p \left( \theta \right|x) \right] =\int q \left( \theta \right) \log \frac{q\left( \theta \right)}{p \left( \theta|x \right)}d\theta$

$=\int q \left( \theta \right) \log \frac{q\left( \theta \right)p \left(x \right)}{p \left( \theta,x \right)} d\theta$

$=\int q \left( \theta \right) \log \frac{q \left( \theta \right)}{p \left( \theta,x \right)} d \theta- \int q \left( \theta \right) \log p \left(x\right) d\theta$

$=\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta-\log p\left(x\right)$

由於對$\log p\left(x\right)$已經使不上力(和$\theta$無關)，所以接下來專心解決$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$，並解我們知道$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$和$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$成正比，所以要盡量把他變小

接下來引入mean field assumption

$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta =\int \prod_i q\left(\theta_i \right) \log \frac{\prod_i q\left(\theta_i \right)}{p\left(\theta,x\right)} d\theta$

$=\int \prod_i q\left(\theta_i \right) \left( \sum_i \log q\left(\theta_i \right) - \log p\left(\theta, x \right ) \right)$

$=\int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \sum_i \log q\left(\theta_i \right) - \log p\left(\theta, x \right ) \right) d\theta$

$=\int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \log q\left(\theta_j \right) - \log p\left(\theta, x \right ) \right) d\theta$ $- \int q\left(\theta_j \right) \prod_{i\neq j} q\left(\theta_i \right) \left( \sum_{i\neq j} \log q\left(\theta_i \right) \right) d\theta$

$= \int q\left(\theta_j \right) \left( \log q\left(\theta_j \right) - \int \prod_{i\neq j} q\left(\theta_i \right) \log p\left(\theta, x \right ) d\theta_{i\neq j} \right ) d\theta_j$$+c$

$= \int q\left(\theta_j \right) \left( \log q\left(\theta_j \right) - E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right) d\theta_j +c$

$= \int q\left(\theta_j \right) \log \frac{\log q\left(\theta_j \right)}{\exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )} d\theta_j +c$

$=KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ] + c$


$KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ]$最小就是讓$\int q\left(\theta\right)\log\frac{q\left(\theta\right)}{p\left(\theta,x\right)} d\theta$最小，也就是讓$KL\left [ q\left(\theta\right) , p\left( \theta \right | x)\right]$最小，目標就完成啦!!!!!!! 哈哈哈哈~~~~~(已瘋)

而可以把$KL\left[ q\left(\theta_j \right), \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right ) \right ]$變得最小的$q\left(\theta_j \right)$不就是$\exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )$本人嗎?

真是可喜可賀!!!

因此Variational Bayes演算法就大功告成了

$q^{*}\left(\theta_j \right) = \frac{1}{Z} \exp \left( E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] \right )$

或者寫成

$\log q^{*}=E_{q \left( \theta_{i \neq j} \right )} \left[ \log(\theta, x)) \right ] - \log Z$

逐個迭代$q \left( \theta_{j} \right )$直到收斂為止

[SAS] 賽仕冷知識(1) missing value 也有大小之分!?

其實平常不太用SAS這個我認為會越寫越笨的東西

題外話，他們倆是關係企業嗎？

但是今天不小心發現了一個它的小秘密!!就是missing value不只一種，然後不同種之間還有大小之分太精細了RRRRRRRRRR~~~

一般來說在SAS裡missing value是用「.」表示，排序上會比有的數字小，所以如果執行 proc sort 程序後，missing value都會在最上面。

除了一般的missing value之外，還有其他種類的missing value。像是「.A」~「.Z」，或是「._」這類特殊的missing value。

根據 SAS的文件，在排序上 「._」是最小的，其次是一般的missing value「.」，再來是「.A」~「.Z」。當然，數字永遠比missing value大。

最後來驗證一下

首先製造一筆資料

data aa; 
input name $ score; 
cards; 
A 60 
B 67 
C . 
D ._ 
E .A 
F .Z
G -10
;

排序一下

proc sort data=aa; 
by score; 
run;

然後來看結果

proc print data = aa;
run;

參考資料
http://support.sas.com/documentation/cdl/en/lrcon/62955/HTML/default/viewer.htm#a000989180.htm#a001221306

[JAVA] Big Data(4) Hadoop Multiple Input

Multiple Input

在map-reduce時，若不同的資料來源要塞給不同的mapper，最後再一起塞進reducer運算，就需要使用Multiple Input 的功能。如下圖，有三個不同的資料來源，先分別進入不同的mapper，然後最後要進到同一個reducer。

首先要先import需要的的類別
org.apache.hadoop.mapreduce.lib.input.MultipleInputs

然後在主程式(main)中寫入下面這一行，告訴電腦你要把哪一筆資料送進哪一個Mapper class的map函數

MultipleInputs.addInputPath(Job名稱, 輸入資料的位址, 格式, Mapper class的名字);

接下來看一個簡單到近乎無腦的例子
假設手上有三筆資料，都包含學號、科目和成績，但是長相就是不太一樣，現在我們要計算各科的平均分數

第一筆

$HADOOP/hadoop fs -cat input_1 | head -3
WUJ-360100;math;56
WPY-802007;math;98
FKT-670008;science;67

第二筆

$HADOOP/hadoop fs -cat input_2 | head -3
{Number:FJB-004150, Subject:math, Score:96}
{Number:QDG-300700, Subject:chinese, Score:90}
{Number:JVY-030140, Subject:chinese, Score:71}

第三筆

$HADOOP/hadoop fs -cat input_3 | head -3
[Number=>ITM-501806; Subject=>science; Score=>82]
[Number=>QBE-003981; Subject=>math; Score=>85]
[Number=>EUJ-017009; Subject=>chinese; Score=>63]

以上三種長相的資料要分別送給三種不同的Mapper中處理，產生(subject, score)的pair然後統一送進一個Reducer做平均數的計算，所以要準備三種Mapper

public static class Map1 extends Mapper
{
      public void map(LongWritable key, Text value, Context con) 
        throws IOException, InterruptedException
      {
              // get the student number
              String stNum = value.toString().split(";")[1];

              // get score
              int score = Integer.parseInt(value.toString().split(";")[2]);
              con.write(new Text(stNum), new IntWritable(score));
      }
}
public static class Map2 extends Mapper
{
      public void map(LongWritable key, Text value, Context con) 
        throws IOException, InterruptedException
      {       
              // "shave" the input value
              String line = value.toString().replaceAll("}", "");

              if(line.contains(",")){
                      // get the student number
                      String stNum = line.split(",")[1].split(":")[1];

                      // get score
                      int score = Integer.parseInt(line.split(",")[2].split(":")[1]);
                      con.write(new Text(stNum), new IntWritable(score));
              }
      }
}
public static class Map3 extends Mapper
{
      public void map(LongWritable key, Text value, Context con) 
        throws IOException, InterruptedException
      {
              // "shave" the input value
              String line=value.toString().replaceAll("[]\\[]", "");


              if(line.contains(";")){
                // get the student number
                String stNum = line.split(";")[1].split("=>")[1];

                // get score
                int score = Integer.parseInt(line.split(";")[2].split("=>")[1]);
                con.write(new Text(stNum), new IntWritable(score));
              }
      }
}

Reducer其實就只需要一個就可以了

public static class Red extends Reducer
{
     public void reduce(Text stNum, Iterable scores, Context con)
      throws IOException , InterruptedException
      {
              int numerator = 0;
              int denominator = 0;
              for (IntWritable v : scores){
                  numerator += v.get();
                  denominator ++;
              }
              int avg = numerator/denominator;
              con.write(stNum, new IntWritable(avg));
      }
}

然後是比較麻煩的主程式

public static void main(String[] args) throws Exception
{
      Configuration conf=new Configuration();
      String[] files=new GenericOptionsParser(conf,args).getRemainingArgs();
      Path inPath1=new Path(files[0]);
      Path inPath2=new Path(files[1]);
      Path inPath3=new Path(files[2]);
      Path outPath=new Path(files[3]);
      FileSystem hdfs = outPath.getFileSystem(conf);
      if (hdfs.exists(outPath)){
        hdfs.delete(outPath, true);
      };

      Job exampleJob = new Job(conf,"example");
      exampleJob.setJarByClass(MpInputExp.class);
      exampleJob.setMapperClass(Map1.class);
      exampleJob.setMapperClass(Map2.class);
      exampleJob.setMapperClass(Map3.class);
      exampleJob.setReducerClass(Red.class);
      exampleJob.setOutputKeyClass(Text.class);
      exampleJob.setOutputValueClass(IntWritable.class);

      MultipleInputs.addInputPath(exampleJob, inPath1, TextInputFormat.class, Map1.class);
      MultipleInputs.addInputPath(exampleJob, inPath2, TextInputFormat.class, Map2.class);
      MultipleInputs.addInputPath(exampleJob, inPath3, TextInputFormat.class, Map3.class);
      
      FileOutputFormat.setOutputPath(exampleJob, outPath);
      System.exit(exampleJob.waitForCompletion(true) ? 0:1);
}

要注意MultipleInputs.addInputPath有沒有把Input和Mapper配對好

最後來看結果(打包部分省略，可以參考這裡)

$HADOOP/hadoop fs -getmerge output_exp output_exp
cat output_exp

science 68
chinese 70
math    68

送上所有JAVA code結束這惱人的一切

import java.io.IOException;
import org.apache.hadoop.conf.Configuration;
import org.apache.hadoop.fs.FileSystem;
import org.apache.hadoop.fs.Path;
import org.apache.hadoop.io.IntWritable;
import org.apache.hadoop.io.LongWritable;
import org.apache.hadoop.io.Text;
import org.apache.hadoop.mapreduce.Job;
import org.apache.hadoop.mapreduce.Mapper;
import org.apache.hadoop.mapreduce.Reducer;
import org.apache.hadoop.mapreduce.lib.input.MultipleInputs;
import org.apache.hadoop.mapreduce.lib.input.TextInputFormat;
import org.apache.hadoop.mapreduce.lib.output.FileOutputFormat;
import org.apache.hadoop.util.GenericOptionsParser;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;

public class MpInputExp
{
  public static class Map1 extends Mapper
  {
        public void map(LongWritable key, Text value, Context con) 
          throws IOException, InterruptedException
        {
                // get the student number
                String stNum = value.toString().split(";")[1];

                // get score
                int score = Integer.parseInt(value.toString().split(";")[2]);
                con.write(new Text(stNum), new IntWritable(score));
        }
  }
  public static class Map2 extends Mapper
  {
        public void map(LongWritable key, Text value, Context con) 
          throws IOException, InterruptedException
        {       
                // "shave" the input value
                String line = value.toString().replaceAll("}", "");

                if(line.contains(",")){
                        // get the student number
                        String stNum = line.split(",")[1].split(":")[1];

                        // get score
                        int score = Integer.parseInt(line.split(",")[2].split(":")[1]);
                        con.write(new Text(stNum), new IntWritable(score));
                }
        }
  }
  public static class Map3 extends Mapper
  {
        public void map(LongWritable key, Text value, Context con) 
          throws IOException, InterruptedException
        {
                // "shave" the input value
                String line=value.toString().replaceAll("[]\\[]", "");


                if(line.contains(";")){
                  // get the student number
                  String stNum = line.split(";")[1].split("=>")[1];

                  // get score
                  int score = Integer.parseInt(line.split(";")[2].split("=>")[1]);
                  con.write(new Text(stNum), new IntWritable(score));
                }
        }
  }
  public static class Red extends Reducer
  {
       public void reduce(Text stNum, Iterable scores, Context con)
        throws IOException , InterruptedException
        {
                int numerator = 0;
                int denominator = 0;
                for (IntWritable v : scores){
                    numerator += v.get();
                    denominator ++;
                }
                int avg = numerator/denominator;
                con.write(stNum, new IntWritable(avg));
        }
   }
  public static void main(String[] args) throws Exception
  {
        Configuration conf=new Configuration();
        String[] files=new GenericOptionsParser(conf,args).getRemainingArgs();
        Path inPath1=new Path(files[0]);
        Path inPath2=new Path(files[1]);
        Path inPath3=new Path(files[2]);
        Path outPath=new Path(files[3]);
        FileSystem hdfs = outPath.getFileSystem(conf);
        if (hdfs.exists(outPath)){
          hdfs.delete(outPath, true);
        };

        Job exampleJob = new Job(conf,"example");
        exampleJob.setJarByClass(MpInputExp.class);
        exampleJob.setMapperClass(Map1.class);
        exampleJob.setMapperClass(Map2.class);
        exampleJob.setMapperClass(Map3.class);
        exampleJob.setReducerClass(Red.class);
        exampleJob.setOutputKeyClass(Text.class);
        exampleJob.setOutputValueClass(IntWritable.class);

        MultipleInputs.addInputPath(exampleJob, inPath1, TextInputFormat.class, Map1.class);
        MultipleInputs.addInputPath(exampleJob, inPath2, TextInputFormat.class, Map2.class);
        MultipleInputs.addInputPath(exampleJob, inPath3, TextInputFormat.class, Map2.class);
        
        FileOutputFormat.setOutputPath(exampleJob, outPath);
        System.exit(exampleJob.waitForCompletion(true) ? 0:1);
  }
}

[機器學習] Paper digest(1) On Multilabel Classification and Ranking with Bandit Feedback

這篇是在Journal of machine learning research上看到的
先附上原文位址

這篇文章提出一個機器學習方法在資訊不充足的情況下做出決策

所謂資訊不足指的是，我們所收集到的資料是未經排序的，但是我們要給的決策建議是卻是排序的。舉例而言，收集到的資料只有某些人的某些特質和他們各自看了哪些電影(partial feedback)，卻不知道這些人對這些電影的喜好程度，我們要如何給出一個排序過的推薦清單?

Loss function for partial information

這篇文章的作者發展的loss function：

$l_{a,c} \big(Y_{t}, \widehat{Y_{t}} \big) =a |Y_{t} \backslash \widehat{Y_{t}}|+\big(1-a\big) \sum_{i \in\widehat{Y_{t}} \setminus Y_{t}} c \big(j_{i}, |\widehat{Y_{t}}| \big)$

其中$Y_{t}$是收集到的資料，$\widehat{Y_{t}}$是建議的決策，因此$Y_{t}$是沒有排序的，而$\widehat{Y_{t}}$是有排序的

$|Y_{t}, \widehat{Y_{t}}|$是指有多少$Y_{t}$裡的element沒有出現在$\widehat{Y_{t}}$裡面，概念上有點像是Hamming loss function；而$\sum_{i \in\widehat{Y_{t}} \setminus Y_{t}} c \big(j_{i}, |\widehat{Y_{t}}| \big)$則是考量了排序，$\widehat{Y_{t}}$中，沒有出現在$Y_{t}$的元素要是排得越前面，計分就越重

舉例來說，如果

$Y_{t}=\left \{ 1,3,8 \right \}$

$\widehat{Y_{t}} = \big(4,3,6\big)$

則

$|Y_{t}, \widehat{Y_{t}}| = 2$ ($Y_{t}$中的$1$, $3$不屬於$\widehat{Y_{t}}$)

$\sum_{i \in\widehat{Y_{t}} \setminus Y_{t}} c \big(j_{i}, |\widehat{Y_{t}}| \big) = 3/3 + 1/3$  ($4$ 和 $6$的權重分別為$3/3$和$1/3$)

Multilabel classification

利用線性模型

$P_{t}(y_{1,t}, y_{2,t}, ...y_{k,t})=P_{t}(y_{1,t}, y_{2,t}, ...y_{k,t}|x_{t})$

$P_{t}(y_{i,t}=1) = \frac{g(-u_{i}^{T})}{g(u_{i}^{T}) + g(-u_{i}^{T})}$

最佳化上述的loss function找最好的決策，上述的loss function也可以寫成

$a \sum_{i=1}^{K}y_{i,t}+\big(1-a\big) \sum_{i \in\widehat{Y_{t}}}\big( c (j_{i}, |\widehat{Y_{t}}| ) - \big(\frac{a}{1-a} + c (j_{i}, |\widehat{Y_{t}}| ) \big)y_{i,t}\big)$

由於$\sum_{i=1}^{K}y_{i,t}$和$\widehat{Y_{t}}$基本上是無關的，所以在計算中只要最佳化

$E\left [l_{a,c} \big(Y_{t}, \widehat{Y_{t}} \big)   \right ]=\big(1-a\big) \sum_{i \in\widehat{Y_{t}}}\big( c (j_{i}, |\widehat{Y_{t}}| ) - \big(\frac{a}{1-a} + c (j_{i}, |\widehat{Y_{t}}| ) \big)p_{i,t}\big)$

就可以了

Others

作者也提出了另一個loss function：

$l_{p-rank, t}(Y, f)=\sum_{i,j\in \widehat{Y}_{t}:y_{i} <  y_{j}}\big(\left \{f_{i}(x_{t} ) < f_{j}(x_{t} )  \right \} + \frac{1}{2}\left \{ f_{i}(x_{t} ) = f_{j}(x_{t} ) \right \} \big)+ S_{t}|Y_{t}\backslash \widehat{Y}_{t}|$

$f$是ranking function而$\left \{... \right \}$是indicator function，$S_{t}$則是$\widehat{Y}_{t}$的長度

作為評估排序表現的依據

[JAVA] JDBC 連接 SQL database

JDBC是一個蠻好用的Package，可以從這裡下載

這裡記錄一些最近摸索的基本的用法，以MS-SQL為例

在使用前記得要(搖一搖?)

import java.sql.*;

0. Connection

String conUrl = "jdbc:sqlserver://portNumber:XX;serverName=XX;databaseName=XX;user=XX;password=*****;";
Class.forName("com.microsoft.sqlserver.jdbc.SQLServerDriver");
Connection con = DriverManager.getConnection(conUrl);

1. Create , table, truncate table

try{
        String query="XXX";//可以是create, drop或是truncate table的sql語句
 
        Statement stmt = con.createStatement();
        stmt.executeUpdate(query);
        stmt.close();
        con.close();

        }catch (Exception e) {
            e.printStackTrace();
        }
    }

2. Insert data

如果只是要插入一筆資料的話，可以把上面的

String query="XXX";

換成

String query="insert into xxx...";

也可以參考prepare statement的寫法，好處是可以不用重複寫sql語句，只要把問號(?)的地方取代成想要插入的值

try{
        String sql="INSERT INTO XXX (column1,column2) VALUES (?,?)";//prepare statement
 
        PreparedStatement pstmt = con.prepareStatement(sql);
        pstmt.setString(1, "插入值1");//第一個?要插入的值
        pstmt.setString(2, "插入值2");//第二個?要插入的值
        pstmt.executeUpdate();
        
        pstmt.close();
        con.close();
        }catch (Exception e) {
            e.printStackTrace();
        }
    }

如果要insert的資料有很多筆，其實一筆一筆塞是有點沒效率的，所以如果要塞入很多筆資料的話，可以先寫入Batch再一次執行

try{
        String sql="INSERT INTO XXX (column1,column2) VALUES (?,?)";//prepare statement
 
        PreparedStatement pstmt = con.prepareStatement(sql);

        for(String s : Data){
            pstmt.setString(1, "要插入的資料");//第一個?要插入的值
            pstmt.setString(2, "要插入的資料");//第二個?要插入的值
            pstmt.addBatch();//寫入Batch
        }

        pstmt.executeBatch();//執行Batch

        pstmt.close();
        con.close();
        }catch (Exception e) {
            e.printStackTrace();
        }
    }

3. 讀料取資

讀取SQL database的資料算是蠻常用的

try{
        String sql="SELECT XXX ...";//Query語句
　　　　　Statement stmt = con.createStatement();
 
　　　　　ResultSet rs = stmt.executeQuery(sql);//Query結果存在這裡
        ResultSetMetaData rsmd = rs.getMetaData(); //取得Query資料

        int numColumns = rsmd.getColumnCount();

        while (rs.next()){//while loop 一筆一筆iterate
        　　　for(int i = 1; i < numColumns+1; i ++){
　　　　　　　　　　　System.out.println(rs.getString(rsmd.getColumnName(i)));//印出資料
　　　　　　　　}
        }

        stmt.close();
        con.close();
        
        }catch (Exception e) {
            e.printStackTrace();
        }
    }

以上是用print當作示範，當然也可以把資料讀下來做其搭更複雜的計算，比如說塞進array裡面bla bla bla..

[Python] 小用法備忘

最近很少跟他打交道，有些小東西雖然基本，但是久不用就會忘記阿阿阿阿阿

lambda

有別於function的寫法，lamda用過即丟，用來建立簡單的function

>>> test = lambda x,y,z: x*y*z
>>> test(1,2,3)
6

filter

用filter函數把不符條件的element濾掉

>>> a = [31,545,21,2,0,231,56]
>>> func = lambda x: x>2
>>> filter(func,a)
[31, 545, 21, 231, 56]

map

用map函數計算每個member的結果

>>> a = [31,545,21,2,0,231,56]
>>> func = lambda x: x+2
>>> map(func, a)
[33, 547, 23, 4, 2, 233, 58]

reduce

用reduce計算$f(f...f(f(x_{0},x_{1}),x_{2})..,x_{n})$

>>> a = [1,2,3,4,5]
>>> func = lambda x,y: x*y
>>> reduce(func, a)
120

range & xrange

用range 函數來產生數列

>>> range(10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

xrange函數用在迴圈中加強效能

>>> for i in xrange(10):
 print(i)

 
0
1
2
3
4
5
6
7
8
9

iterate寫法

iterate寫法是python特有的解決方案，和for loop類似

[i*5 for i in xrange(4)]
[0, 5, 10, 15]